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Abstract 

The aim of this review is to develop the kinetic theory of phonons in classical parti¬ 
cle chains to a point which allows comparing the kinetic theory of normally conducting 
chains, with an anharmonic pinning potential, to the kinetic theory of the anomalously 
conducting FPU chains. In addition to reviewing existing results from the literature, we 
present a streamlined derivation of the phonon Boltzmann collision operators using Wick 
polynomials, as well as details about the estimates which are needed to study the effect 
of the collision operator. This includes explicit solutions of the collisional constraints, 
both with and without harmonic pinning. We also recall in detail the derivation of the 
Green-Kubo formula for thermal conductivity in these systems, and the relation between 
entropy and the Boltzmann H-theorem for the phonon Boltzmann equations. The focus 
is in systems which are spatially translation invariant perturbations of thermal equilib¬ 
rium states. We apply the results to obtain detailed predictions from kinetic theory for 
the Green-Kubo correlation functions, and hence the thermal conductivities, of the chain 
with a quartic pinning potential as well as the standard FPU-/3 chain. 

Contents 

1 Kinetic scaling limit for weakly anharmonic chains [2 

1.1 Introduction. [2] 

1.2 Free motion of phonons . 0] 

1.2.1 Effect of boundary conditions in finite non-periodic chains. [ 6 ] 

1.3 Particle chains with anharmonic perturbations . [ 6 ] 

1.4 Choice of initial data. . 

1.5 Green-Kubo formula. M 

2 Phonon Boltzmann equation of spatially homogeneous anharmonic chains 1191 

2.1 Identifying the collision operator. [19] 

2.2 Solution of the collisional constraints. [25| 

2 . 2.1 C 3 = 0 for nearest neighbour dispersion relations. 021 

2 . 2.2 C 4 for nearest neighbour dispersion relations. [27| 

2.2.3 Integration of the collisional constraints in FPU chains. 05] 

2.2.4 Integration of the collisional constraints for onsite nonlinearity .... [30| 


*E-mail: jani . lukkarinenShelsinki .f i 


1 














3 Energy transport in the kinetic theory of phonons 1331 

3.1 Entropy and H-theorem of the phonon Boltzmann equations. [33] 

3.2 Steady states. |35| 

3.3 Green-Kubo formula and the linearized Boltzmann equation . [36| 

3.4 Kinetic theory prediction for thermal conductivity in chains with anharmonic 

pinning . [371 

3.5 Anomalous energy conduction in the kinetic theory of FPU chains. [JT] 

4 Concluding remarks 1441 


1 Kinetic scaling limit for weakly anharmonic chains 


1.1 Introduction 

Kinetic theory describes motion which is transport dominated in the sense that typically the 
solutions to the kinetic equations correspond to constant velocity, i.e., ballistic, motion in¬ 
tercepted by collisions whose frequency is order one on the kinetic space-time scales. The 
constant velocity part of the transport, the free streaming motion, can arise via several differ¬ 
ent microscopic mechanisms, the obvious case being free motion of classical particles. Here, 
we are interested in evolution of energy density in chains composed out of locally interacting 
classical particles, in the limit where the interactions are dominated by a stable harmonic po¬ 
tential. For such models, the energy transport is dominated by the free streaming of phonons 
which correspond to eigenmodes of the motion generated by the harmonic potential. 

Our goal is to derive and study the properties of kinetic theory given by a phonon Boltz¬ 
mann equation. Even with most generous estimates for the terms neglected in its derivation, 
the Boltzmann equation can become an exact description of the energy transport only in a 
kinetic scaling limit which produces a total scale separation between the free streaming mo¬ 
tion of phonons and the collisions induced by the perturbation. For the above particle chains 
this amounts to studying lengths and time intervals both of which are proportional to A“^, 
where A is the strength of the anharmonic perturbation, and then taking the limit A —>• 0. 

In particular in one dimension, i.e., for particle chains, it is difficult to estimate reliably 
the error terms associated with the approximation of the original motion by the phonon Boltz¬ 
mann equation. At present, no such full mathematical analysis has been achieved. However, 
there are instances in which the one-dimensional phonon Boltzmann equation has proven to 


be useful: in Section 3.4, the thermal conductivity predicted by the Boltzmann equation of 
systems with an anharmonic onsite potential is compared to numerical simulations, and the 
results are found to agree within numerical accuracy for sufficiently weak couplings. 

Although just one example, the agreement is significant since it shows that, under proper 
conditions, the kinetic theory from the phonon Boltzmann equation can produce meaningful 
predictions even for one-dimensional transport problems. As will be apparent later, com¬ 
puting these predictions requires some additional effort compared to the standard rarefied 
gas Boltzmann equation. The effort is rewarded in a prediction which has no adjustable 
parameters and hence can be compared directly with experiments and numerical simulations. 

The drawback of the phonon Boltzmann equation is that, a priori, it only describes the 
lowest order effects of the perturbation. Hence, it might not capture all processes relevant to 
the energy transport. In addition, since it is not known which precise conditions guarantee 
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that the corrections are truly vanishing in the kinetic scaling limit, one has to leave the option 
open that even the lowest order effects are not being described accurately by the equation. 

A discussion about the physical and mathematical conditions under which the phonon 
Boltzmann equation should be an accurate approximation is given in [T] , while [ 2 ] could serve 
as a textbook reference on physics of phonons. The present rigorous results on linear and 
nonlinear perturbations of wave-like evolution equations support the basic principle that for 
sufficiently dispersive systems, in particular for crystals in three or higher dimensions, the 
derivation of the phonon Boltzmann equation yields an accurate approximation to the be¬ 
haviour of the system at kinetic time-scales with weak coupling. The argument is particularly 
convincing for spatially homogeneous systems and for linear perturbations such as when the 
particle masses have a weak random disorder laiii, and for the closely related quantum mod¬ 
els, random Schrodinger equation 13 E] and the Anderson model [7] . In addition, it has been 
proven that the Wigner transforms of unperturbed wave-equations do, in great generality, 
converge to the collisionless Boltzmann equation [SIE], although acoustic modes may require 
somewhat more refined treatment |10j . Nonlinear perturbations are more intricate, and we 
refer to nulls] for recent results on the rigorous analysis of their kinetic scaling limits. 

Quite often the addition of the collision operator to the transport equation leads to loss 
of “memory” at each collision and eventual equilibration via diffusion: this would correspond 
to normal energy transport by phonons. However, if the transport is anomalous, it is highly 
likely to be reflected also in the solutions of the corresponding kinetic equation. An example 
of such behaviour is given by the FPU-/3 chain. The aim of this contribution is to develop 
the general kinetic theory of classical particle chains to the point which allows to compare 
the kinetic theory of normally conducting chains with an anharmonic pinning potential and 
the anomalously conducting FPU chains. 

The first section concerns the derivation of the phonon Boltzmann equation in the one¬ 
dimensional case, for general dispersion relations and third and fourth order perturbations of 
the harmonic potential. We begin with the infinite volume harmonic model in Sec. | 1 . 2 [ and 
make a brief comment on the finite volume case and the relevance of boundary conditions in 


Sec. 1.2.1 The anharmonic perturbations are discussed in Sec. 1.3 where we also explain our 


choice for the definition of the related energy density and derive the related energy current 
observables. 

The role of the proper choice of initial data for the derivation of the Boltzmann equation 
cannot be stressed enough. We describe our choices in Sec. 1.4 and introduce the notations 
needed in the analysis of the cumulants of the phonon fields. We conclude the first section 
with a derivation of the Green-Kubo formula for the particle chains, under the assumption 
that energy density satisfies a closed diffusion equation. (This would correspond to a case 
with only one locally conserved field.) 

The first section concerns fairly general particle chains, but in order to derive the phonon 
Boltzmann equation, further restrictions are required. In Sec. 2.1, we do the derivation in 


detail for the FPU chains, and sketch it for the onsite anharmonic perturbations. The rest of 
the section concerns only chains with stable nearest neighbour harmonic couplings. We explain 
why the nearest neighbour couplings do not produce any collisions at the kinetic time-scale 
from the third order terms of the potential (Sec. 2.2.1), and show that only phonon number 
conserving collisions appear from the fourth order terms (Sec. 2.2.2). To use the collision 
operator, one needs to resolve the related collisional constraints. This is done for FPU-chains 


in Sec. 2.2.3 and for the chains with a harmonic pinning potential in Sec. 2.2.4 


The final section is devoted to the analysis of the resulting class of phonon Boltzmann 
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equations. We first show how an assumption of an increase of microscopic entropy in homo¬ 
geneous models leads to an H-theorem for the phonon Boltzmann equation. In Sec. 3.2 it 
is explained how the H-theorem can be used to classify all steady states of the Boltzmann 
equation. The kinetic description derived in Sec. 1.5 for the time correlation function in the 
Green-Kubo formula is applied to these models in Sec. |3.3[ This results in explicit predic¬ 


tions for the leading asymptotics of the decay of the correlations. The decay is found to be 
integrable for the onsite perturbation, and we make a comparison of the resulting thermal 


conductivities in Sec. |3.4[ We discuss the corresponding prediction in FPU models in Sec. 3.5 
The result is decay which predicts anomalous conduction in the FPU chains. 

More comments and discussion about the results are given in Sec. before the Acknowl¬ 
edgements and References. 


1.2 Free motion of phonons 

Phonons correspond to eigenmodes of harmonic Bravais lattices. The general case includes 
multi-component phonon fields which may arise both from the dimensionality of the lattice— 
the evolution equation could describe displacements from some reference lattice positions and 
thus have d components in d dimensions—or from the reference cell of periodicity containing 
more than one particle. An example of the latter is provided by the one-dimensional chain with 
alternating masses: then the reference system remains invariant only under shifts over two 
lattices sites which corresponds to a Bravais lattice with two particles per cell of periodicity. 
For notational simplicity, let us here only consider the case of one-component classical phonon 
fields and assume that all particles have a unit mass (this is always possible to achieve for 
one-component fields by choosing a suitable time-scale). 

Then the state of the system at time t is determined by (qxit),Pxit)), x £ Z, which satisfy 
the evolution equations 


qx{t) =Px{t), Px{t) =- y)qy{t). (1) 

y 

The equations are of a Hamiltonian form, corresponding to the Hamiltonian function 

H{p, q) = '^\pl+ \qxa{x - y)qy . (2) 

x£lj x^y&, 

We suppose that the harmonic interactions have a short range. For instance, only finitely 
many a{z) are different from zero, or |a(^)| decreases exponentially fast to zero as \z\ —>• oo. 

The evolution equation Q can be solved by using Fourier-transform. We define the 
Fourier-fields as q{t,k) = Here k £ T 

where T denotes the one-torus which we parametrize using the interval [—|, and then 
identify its endpoints, — ^ and The evolution equation for the Fourier-fields then reads, in 
a matrix form. 


±(q{t,k)\^( 0 1\ fq{t,k)\ 

dt\p{t,k)J \-a{k) oj \pit,k)J 

For each k, the equation is easily solved by Jordan decomposition of the matrix on the right 
hand side. This provides two independently evolving eigenmodes at{k,a), cr = ±1, as long 
as d{k) / 0. For those k with d{k) = 0, the solution is given by constant speed increase: 
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p{t, k) = p(0, k) and q{t, k) = g(0, k) + tp{0, k). Typically, this can happen at most at a finite 
number of points whose neighbourhoods might require special treatment in the kinetic theory 
of phonons. 

The eigenmode fields at{k,a), which we call phonon modes of the harmonic chain, satisfy 
the evolution equation 

^at{k, a) = -iauj{k)at{k, a) , (4) 

at 

where ui{k) = ^ya{k). Thus the square root of the Fourier-transform of a determines the 
dispersion relation oj of the phonons, and we arbitrarily choose the principal branch of the 
square root so that always Rea;(A:) > 0. There are several possible ways to normalize the 
eigenmodes. Here we employ the following choice, standard to phonon physics, 

atik,a)= , {uj{k)q{t,k)+ iap{t,k)), a£{±l},k€T. (5) 

^/2u}{k) 

It is straightforward to check that these fields indeed satisfy the evolution equation Q 
whenever u}{k) ^ 0. In addition, if we choose a solution to @ and then define 

q{t,k) =—^='^at{k,a) and p{t, k) = '^{-ia)at{k, a ), (6) 

then {q{t, k),p{t, k)) yields a solution to 0. Given some initial field oq, the solution of Q is 
straightforward, yielding at{k,a) = u). Therefore, if we let qq to be determined 

by the initial data (q(0, k),p(0, k)), we can now obtain the solution to the original evolution 
problem ^ for every k for which a{k) ^ 0, just by inserting at{k,a) = e~^^‘^^^^^aQ(k,a) in 

(©. 

For later use, it is also important to assume that the harmonic interactions are stable, 
i.e., they do not have any solutions which increase exponentially in time. By the above 
discussion, this is equivalent to assuming that the Fourier-transform of a is pointwise non¬ 
negative, a{k) > 0 for all k € T. A standard example is given by the nearest neighbour 
interactions for which a(z) = 0 if \z\ > 1. A short computation reveals that the only stable 
nearest neighbour dispersion relations are given by ijj{k) = a;(l — 2(Icos(27r(fe—fco)))^/^ for 
some tJ > 0, |(5| < ^ and ko G T. Since for phonons we also require that a{z) are real, we can 
also always choose /cq = 0 above. 

The parameter w determines the average angular frequency of the oscillations, the maxi¬ 
mum value begin -|- 2|(I| and the minimum — 2|(I|. If |5| < the function u}{k) is 
a thus an even, real analytic, function. These systems describe optical phonon modes, and 
one can also think of the system as having harmonic pinning. If <5 < 0, it is possible to 
make a change of variables which reverses its sign to positive: this is achieved with a shift 
by ^ in the Fourier variable, in other words, by using {{—l)^qx, i—l)^Px) as the new lattice 
fields. Enforcing this convention has the benefit that then the slowest evolving fields are also 
spatially slowly varying, corresponding to small values of k. Therefore, we only consider the 
cases ko = 0 and 0 < d < | in the following, excluding also the degenerate constant dispersion 
relation corresponding to d = 0. 

If (5 = ^ one has uj{k) = aJ\/2| sin(7rA:)|, with an | A:[-singularity near origin. These systems 
describe acoustic phonon modes and the seemingly innocuous singularity has important im¬ 
plications for the energy transport properties of the system. The FPU-chains belong to this 
category. 
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Finally, let us point out that in the stable case, Ljj{k) > 0, the phonon fields defined by Q 
are directly connected to the energy of the harmonic system by ^ f dkLo(k)lat(k, = 
H{p{t), q{t)). Hence, in analogy to quantum mechanics, the complex phonon fields at{k,a) 
can be thought of as wave functions of phonons each of which carries an energy uj{k). 

1.2.1 Effect of boundary conditions in finite non-periodic chains 

For notational simplicity, the above discussion uses formally an infinite lattice setup. Its 
proper interpretation is to consider the infinite system to be an approximation for a finite, but 
large system. This identification has an important mathematical consequence: the standard 
methods of £^-spaces do not immediately apply to the infinite system. For instance, typical 
samples from infinite volume Gibbs states are only logarithmically bounded at infinity m- 
Hence, should the need arise, it is better to go back to the finite systems to resolve any 
possible issues about existence of solutions or to handle singularities. (For discussion about 
the existence and properties of the infinite volume dynamics, see Haim.) 

In a finite system, the choice of boundary conditions begins to play a role. The above 
computations involving Fourier-transforms can be given an exact correspondence by choosing 
a box with periodic boundary conditions. In this case, the wave evolution can move energy 
around the lattice unimpeded, as indicated by the solutions which are given by multiplication 
operators in the Fourier-space. One can find more details about the correspondence, for 
instance, in | 12 ) . 

If the finite lattice does not have periodic boundary conditions, the boundaries will influ¬ 
ence the transport of phonons. For boundary conditions which preserve energy, such as Dirich- 
let and Neumann, the effect can be understood as a reflection of phonons at the boundary. 
It is also possible to combine these with partial periodic transmission through the boundary. 

The effect of standard heat baths, such as Langevin heat baths, on the kinetic equations is 
more difficult to analyse. These will introduce also absorption and source terms to the phonon 
evolution equations, but since the heat baths are typically acting only on a few particles of the 
chain, their coupling to the wave evolution is nontrivial and somewhat singular. (In general, 
it is difficult to affect the long wavelength part of the evolution by any local change to the 
system.) Unfortunately, a systematic analysis of these effects appears to be missing from the 
current literature. 

1.3 Particle chains with anharmonic perturbations 

For the discussion about the effect of local nonlinear perturbations to the above wave equation, 
let us introduce, in addition to the harmonic coupling function a : Z —)• M, its second and 
third polynomial counterparts, 03 : Z^ —>■ M and 04 : Z^ —)• M. These are assumed to be 
“small” in the sense that for large microscopic times the evolution is dominated by the linear 
term. Explicitly, we take the perturbed evolution equations to be 

qx{t) =Px{t), 

Px{t) = - '^a{y)qx-y{t) - ^ a3{yi,y2)qx-yi{t)qx-y2{t) 

y yi,y 2 

- a4iyi,y2,y3)qx-yAthx-y2 {t)Qx-y3 it) • (7) 

yi,y2,y3 
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This corresponds to Hamiltonian evolution with the interaction potential 

r( 9 ) = ^E a{y)qxqx-y + ^ ^ a3{yi,y2)qxqx-yi{t)qx-y2{t) 

x,yi,y2 

+ ^ E ®4(yi) y2i yd,)qxqx—y\it)qx—y2it)qx—yo,(t ), 

x,yi,y2,y3 

supposing, as we will do here, that the following symmetries hold: 
a{y) = a{-y ), 

a3{yi,y2) = a3{y2,yi) = a3{-yi,y2 - yi), 

ai{yi,y2,y3) = a4(y2,yi,y3) = a4(yi,y3,y2) = a4(-yi,i/2 - yi,y3 - yi) • 


(8) 


(9) 

( 10 ) 

( 11 ) 


The symmetries arise from assuming that the coefficients correspond to a generic label- 
exchange symmetric, translation invariant potential. We have also used here the fact that any 
permutation can be expressed as a product of adjacent transpositions. Therefore, it suffices 
to check symmetry with respect to the transpositions to get full permutation invariance. 

For example, the standard FPU-chain is given by 

Vppviq) = ^U{qx-qx-i), U{r) = + X 3 ^r^ + , (12) 

X 

and it satisfies the evolution equation qx = U'{qx+i — qx) — U'{qx — qx-i)- A brief computa¬ 
tion, using the periodicity, shows that the FPU chain corresponds to choosing the coefficient 
functions above so that their nonzero values are 


a{y) 


a3{yi,y2) 


a4(y) 


l_ 2 l 2 , 2 / = 0 

-uj < , 

2 1-1, ye{±l} 




2 , 

1 , 

- 1 , 


yi = y 2 = 1, yi = -i, 2/2 = 0 , or yi = 0 , 2/2 = -i 
2/1 = 2/2 = -1, 2/1 = 1 , 2 / 2 = 0 , or 2/1 = 0, 2/2 = 1 

2/1 = 2/2 = 2/3 = 0 

±yG {(0,1,1), (1,0,1), (1,1,0)}, 

±yG 1(1,1,1),(1,0,0),(0,1,0),(0,0,1)} 


(13) 

(14) 


(15) 


For later use, let us also record their Fourier-transforms. Using the shorthand notations 
Pi = 27rki, i = 1,2,3, they are given by 


a{ki) = a;^(l — cospi) = 2uj^ sin 


23 (^ 1 , ^ 2 ) = iA 32 '* sin 

(Xi{ki,k2,k3) = A 42 Re|^]^(l -e 


1—2 „;„2 Pf 
2 ’ 

3 . Pi+P 2 . Pi . P 2 
sm — sm — , 
2 2 


-\pi^ 


1 o 4 ■ PI+P 2 +P 3 
= —A 42 sm--- 


n- 


Pe 
sm — . 


e=i i=i 

Onsite perturbations also belong to the above category. For instance, 

Vos(g) = ^ yxa{x-y)qy + ^V{qx), V {q) = X3^q^ + X4^q^ , 

x,y&Z X 


(16) 

(17) 

(18) 


(19) 
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has 


a3(?/i,?/2) = A3l(yi = y2 = 0), (20) 

a4(2/i,2/2,y3) = XiMvi = y2 = 2/3 = 0). (21) 


Both Fourier-transforms are thus constant: a^{ki,k 2 ) = A 3 and 04 (^ 1 , ^ 2 , ^ 3 ) = A 4 . 

After these examples, let us come back to the general case, equations Q. Taking a 
Fourier-transform yields 


^g(t, fe) =p{t,k ), 

^p(t, k) = -a{k)q{t, k) - [ dk[dk 2 djik - k[ - k 2 )d 3 {k[,k 2 )q{t, k[)q{t, fe^) 
di Jt2 


/T3 


3 

dk[dk 2 dk' 2 , 6j{k — A:()a 4 (k') q{t, k'l ), 

£=l 


£=1 


( 22 ) 


where we have denoted the periodic 5-function by 5 t. We will drop the subscript in the 
following and merely use 6. This implies that the phonon fields, defined still by (§ , now 
satisfy the evolution equation 


—at(A:o,ao) = -iaooj{ko)at{ko,ao) 


io-Q / dkidk2S{ko - ki - k2)^3{ko,ki,k2)at{ki,ai)atik2,a2) 




/T 2 


„ 3 3 

i(To E L dkidk 2 dk 36 {ko — E ke)^4{ko,ki,k2,k3)Y\at{ke,(Ji) 

Cl,a'2,0'3£{±l} ^=1 ^=1 


(23) 


where the interaction amplitude functions are 


^>3(A:o,fei,/i:2) = S3 (A:i,A:2)TT- - -(24) 

£=o (2w(fe£)) 2 

3 1 

$ 4 (^ 0 , ki,k2, k3) = d4{ki, k2, k3) TT-^ . (25) 

e=o 2 

Since the 5-function enforces the sum of the integration variables to be equal to ko, we can 
thus employ the following definitions for the FPU-chains 


2 


$3(/co,A;i, A: 2 ) = iA32Jcj ^Y{gike), 

e=o 

(26) 

3 

^>4(fco, ki, k2, k3) = -2A4w“^ q giki) , 

(27) 


e=o 


where g{k) = sign(fc)y| sin(7rA:)| for \k\ < 

To study energy transport, it is necessary to split the total energy into local components. 
The harmonic energy, corresponding to the energy of free phonons, is most conveniently 
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distributed using a form which is symmetric in Fourier-space. By the definition Q, we have 
for any ki , ^2 


u:{ki)uj{k 2 ) ^ a(A:i, -a)a{k 2 , a) = u{ki)ui{k 2 )q{ki)q{k 2 ) + p{ki)p{k 2 ) ■ (28) 

cr=dil 

Thus, if we denote the inverse Fourier transform of w by oj, we can define 

Hf\p,q) = ]:pl + ]:{^^{y)qx-^ = I Y. [ dke^^^^^'‘^+'^^'><^2{k)a{ki,-a)a{k2,a), 

(29) 


where ^ 2 {ki-,k 2 ) = ^Juj{ki)uj{k 2 ). Since u is real and symmetric, we have here uj{x) G 
M. Thus > 0 and Ylx 9 ) i® equal to the harmonic part of the total energy, 

i YlxPx + 5 Ylxy — y)qy Thus this choice allows distributing the positive total 

harmonic energy into positive contributions localized at each lattice site. 

For the anharmonic terms of the potential energy we use the form suggested by the 
notation in ([^, and define the local energy at site x by 


Hx{p,q) = H^^\p,q) + ^ E c^3(yi) y2)qxqx—y\qx—y2 

yi,y2 

+ i E aA{yi,y2,y3)qxqx -yiQx-y2qx-y3 

yi,y2,y3 

= - V / dA;e‘2^^(^i+^2)$2(A;)a(A:i,-fT)a(A:2,f7) 

1 r ^ 

betSp-'T’ I.\ 

1 r ^ 

+ 4 E 


(30) 


f7e{±i}4 


Then clearly Hx{p, <?) = 5 YlxPx + ^( 9 ) = H{p, q), and and is a local function at x (it 
depends mainly on fields near the point x). In addition, Hx is translated just like the fields 
if these are shifted by xq: if we choose some xq G Z and define qx = qx+xQ and p^ = Px+xo-, 
then Hx{q,p) = Hx+xo{q,p)- 

There are various possible choices for how to split the total energy into local energy density. 
The above definition of energy density is perhaps not the most standard one, but it has two 
appealing features for those systems which are dominated by phonon transport. First, it has 
simple algebraic dependence on the phonon eigenmode fields. Secondly, its dependence on 
the position variable x is located entirely in the Fourier-factor. This allows a definition of 
a current observable with a simple dependence on the particle interactions and the simplest 
choice for the associate discrete derivative. 

Explicitly, consider some k ^ Q and any xi, X2 G Z with xi < X 2 - Setting y = X 2 — xi, we 
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then have y > 0 and 


y 

i27ra:/c _ i27ra:i/c 




X=Xl 


= e—“ ^ = e 

y'=o 


k _ i27ra:ifc 


I — gi27r(y+l)fc 


1 - e 


'l27Tk 


2 sin( 7 rA:) L 


gi27ra:i/c _ gi27r(a:2 + l)A; 


( 31 ) 


Hence, if we define Jx-i,x{t) by replacing in (30) every factor (t^)” , where 

k = by {W(,(^t{k(.,(T()y\ we have for all xi,X 2 € Z with xi < X 2 that 


dt ( Hx{q{t),p{t)) I = - Jx 2 ,x 2 +i{t ). 


(32) 




Since xi,X 2 are arbitrary, we can interpret Jx,x+i{t) as the energy flux from site x to x + 1 
at time t. 

To be precise, the above construction can only be used inside the integral if 7 ^ 0. To see 
why the resulting current observable still works, let us for the moment regularize the system 
by replacing the infinite lattice Z by a finite periodic lattice of length L 3> 1. Then ki G Z/L, 
and either k = G Z or the distance of k from Z is at least 1 /L. Let us thus separate all 

terms with fe G Z in the definition of Hx, and use the periodic identification. Then the terms 
correspond to fc = 0. Since for every x on the periodic lattice e'^™^l(fe=0) = ^ 
these terms sum to H{p,q)/L, i.e., to the total energy density. By conservation of energy, 
these terms do not contribute to dtHx{q{t),p{t)). The remaining terms converge to the 
claimed result as L —)• 00 , provided that the integrals over ¥"■ are understood as principal 
value integrals around the subset with k = 0. 

In summary, the definition 
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Jx-l,x = ^Y / ^k^2{k) ^ 

2 Jf2 2 sm(7rK) 


1Q ITT k _ 

^2'Kxk 


dt {at{ki, -a)at{k2,cr)) 


k=k\+k2 


13 -Its 2 siniTTk) 


■,p-'mk 

^2'Kxk 


7e{±i}3 


dt '[\.at{ki,ae) 


k—Yi—1 1 


1 


;„-i7rfc 

^i2TTxk 


+ - y / dk<^4{k) , 

2sm(7rA;) 


dt '[\.at{ki,ae) 


(33) 




yields a current observable associated to the energy density (30) which satisfies the lattice 


continuity equation (32). Let us point out that no spurious imaginary part is created by this 


choice: since at{k, a)* = at(—k, —cr), $„(A:)* = <!>„(—/c), for all n, and uj{—k) = uj{k), we have 
always Jx,x+i{t)* = Jx,x+i{t) and hence the current observable is real-valued. We may also 


simplify the prefactor above by using 


= I (1 + icot(7rA:)). 


2 sin(7rAi) 

The above current observable also behaves as expected with respect to spatial translations 
of the lattice. If we choose some xq G Z and define qx = qx+xo and px = Px+xo-, then the a-field 
of the translated configuration {q,p) satisfies dt{k, a) = e^‘^'^^°^at{k, a). Thus each of the three 
terms in the sum defining the current in (33) will acquire a factor with k = X)£=i ki- 

Hence, the current of the translated configuration satisfies Jx,x+i = Jx+xq,x+xo+1: similarly 
to what was proven earlier for the energy density observables Hx- 
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1.4 Choice of initial data 

We consider the above evolution equation with random initial data. The distribution of 
initial data plays as important a role as the choice of the coupling functions. For instance, 
there always are degenerate initial data which do not thermalize. Assume, for instance, that 
g(0) is chosen as any local minimum or maximum configuration of the potential V and let 
the particles be initially at rest, p(0) = 0. Since then VV((7(0)) = 0 this conhguration is a 
stationary solution of the evolution equations. However, it need not be stable, and it can well 
happen that even a minute perturbation of the initial data will take the asymptotic state of 
the system far from the original stationary state. 

On the other hand, for normally conducting systems sufficiently chaotic initial data should 
lead to thermalization of the state. The precise mechanism and classihcation of what kind 
of initial data should have this property, and in which precise mathematical sense, is a long¬ 
standing open problem in mathematical physics. The most difficult part of the problem seems 
to be to control the beginning of the thermalization process, the microscopic thermalization, 
where the state of the system approaches local stationarity. 

However, many important physical properties of the system do not require full understand¬ 
ing of the thermalization process. For instance, for systems with normal heat conduction, the 
Green-Kubo formula allows to compute the thermal conductivity of the system. As shown in 
the next subsection, the formula involves studying the time-correlations of the system start¬ 
ing it from a thermal equilibrium state or, equivalently, to study the evolution of correlation 
functions for small perturbations of the equilibrium state. Hence, thermal equilibrium states 
and their perturbations form an important class of initial data, and our assumptions should 
allow for at least these. 

To arrive at the kinetic theory of phonons, we need to make some rather specific assump¬ 
tions about the initial data. The principle here is that the states are chosen to mimic states 
produced by typical microscopic time-averages. Note that this concept involves a choice of 
scale, the time-scale over which the time-averages are taken. These questions touch upon the 
tough question of microscopic thermalization, and we do not wish to speculate further about 
them here. Instead, let us postulate a number of assumptions on the initial data based on 
reasonable assumptions about the physical characteristics of the time evolution. 

We assume here that the initial data for the system has the following properties. 

1. Random: We suppose that the initial data (9a;(0),Pa;(0)) is randomly distributed. The 
initial randomness makes also the conhgurations {qx{t),Px{t)) at later times random 
variables. 

2. Chaotic: We suppose that initial distribution of particles in regions which are far apart 
are almost independent. More precisely, we suppose that all correlation functions of the 
initial data decay fast (on the microscopic scale). The speed is assumed to be at least 
so fast that the correlation functions are .^2-summable. 

3. Spatially homogenized: We assume that there is a scale which is large in micro¬ 
scopic units, such that the correlation functions are nearly invariant under translations 
of lengths less than e~^. 

The randomness in the initial data makes also the phonon helds at{k, a) random variables. 
In principle, the field is not dehned for k such that uj{k) = 0, but we assume these to have been 
suitably regularized. Consider for instance the case in which the interactions depend only on 
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the differences qx — Qx-i, such as for the FPU-beta chains. Given any initial data (on a finite 

periodic chain), we can regularize it before defining the phonon fields without altering the 

evolution equations in any way. We move from {q,p) to {q,p) by first setting p = ^ 

and q = ^ Yhx then we define qx{t) = qx{t) — q — tp and Px{t) = Px{t) — P- Under 

these assumptions, the total momentum is conserved, and thus YlxP^i^) = 0 = 

all t. In addition, {qx{t),Px{'t)) sue a solution to the original evolution equations. We take 

at{k, a) to be defined via these regularized fields which guarantees that 0^(0, cr) = 0 for all t. 

As we are interested in the evolution of the energy density (30), it suffices to study the 
correlation functions of at{k,cr). For instance, to study the energy density averaged over the 
initial data, one needs the correlation functions up to order 4. 

To control the correlation functions of the above type of chaotic, spatially homogenized 
states, it is often better to consider the cumulants of the fields rather than their moments. 
For instance, consider the field px at two points xi,X 2 which are far apart. By the chaoticity 
assumption, then ei = and 62 = Px^ should be nearly independent random variables, and 
hence (e”e™) ~ (ei)(e™), for n,m ^ 0. Unless one of the particles is frozen to its initial 
position, these moments are not zero. However, the corresponding cumulant is then nearly 
zero whenever both n,m ^ 0. 

Hence, for two asymptotically independent regions, it is the cumulants, not moments, 
which will vanish in the limit of taking the regions infinitely far from each other. We consider 
here only initial states for which this decay is absolutely square summable in space, “£2- 
clustering” for short. This particular assumption allows for an easy identification of the 
(distribution type) singularities of the Fourier-transform for any spatially homogeneous initial 
data. 

To make the computations manageable, we rely here on the notations and basic re¬ 
sults used in m- In particular, for a random field x G Z, and any sequence J of 

n lattice points the shorthand notations {il^{x)'^), K['ip{xj)], and :'ip{x)'^: refer to the expec¬ 
tation (nr=i V’(3^j ))) to the cumulant «;[^/^(xjJ,..., V’(xj„)], and to the Wick polynomial 
:nr=iV’(xj. ):, respectively. In addition, we define 'ip{x)‘^ = I if J is an empty sequence. 
Similar notations will be used for random fields on other label sets in addition to the lat¬ 
tice points; for instance, for the fields a{k,a) and any sequence J = {ki,cri)^^^ we write 
K[a{k,a)j] = K[a{ki,ai),... ,a{kn,(Tn)]- 

We employ here the following basic relations between the above constructions. The proofs 
of these results can be found in many sources, in particular, in |15j using the above notations. 
The moments-to-cumulants formula states that for any sequence J 


{i’ixY) = X] n Y^’Ya)] , (34) 

Trer{.J) A&n 

where V{J) denotes the collection of partitions of the sequence ■'0 For a partition vr G 
'P(I), let us call the subsets A G vr clusters or blocks. We also recall that the cumulants 
are permutation invariant and multilinear, which means that they are linear in each of the 
arguments separately. Thus cumulants of Fourier-transforms can be easily expressed as a 
linear combination of the cumulants of the original field. 

^Some care is needed in the interpretation of the moments-to-cumulants formula in order to get all combi¬ 
natorial factors correctly. It is safe to think that the elements in the sequence J of length n are labelled by the 
set 7„ = {1,2,. ..n}, and the collection of partitions V{J) refers to the collection of standard set partitions 
of In- Therefore, if 7 £ A £ tt £ P{J), we have I £ In, xi denotes the value of x at the 7:th position of the 
sequence J, and xa denotes the subsequence {xe)e£A of J. 
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Wick polynomials of random variables allow simplifying the moments-to-cumulants ex¬ 
pansions by removing all those partitions vr from the sum which contain a cluster internal 
to the Wick polynomial part of the expectation. Explicitly, consider some L > 1 and a col¬ 
lection of L -|- 1 index sequences J' and for £ = 1,..., L. Then, for the merged sequence 
I = J' + have 

(^ip{xy'Yl:i^{xy^:\= Yl{K['il^{xA)]HA(;t Je\/l)) . (35) 

' i=i ' w&ni) agtt 

In particular, the formula implies that cumulants and Wick polynomials satisfy the following 
relation: If the sequence J is non-empty and has xi as its first element, then 

yi’ixj)] = {i’ixi) , (36) 


where J \xi denotes the sequence which is obtained when xi is removed from J. By per¬ 
mutation invariance of the cumulants, the formula holds also if xi is replaced by any other 
element in J. 

Let us also point out that, since the cumulants are multilinear, their time-evolution can 
be written down fairly easily for the present kind of hybrid dynamics, where all randomness 
is contained in the initial data. Namely, by using a straightforward “telescoping” argument 


and (36), it is clear that for any sequence J 


dtK.[il)t{xj)] = YidtMxe) ■■'ipixy'^^^:) 

t&J 


(37) 


For the manipulation of terms such as dt'4^t here, it is possible to move back and forth between 
standard and Wick polynomials by using the following identities: 


ucJ 


(38) 


and 


■i’ixy: = Y 11 , ( 39 ) 

U<ZJ TT&P{J\U) A&TT 

where |7r| denotes the number of clusters in the partition vr. 

Let us for the moment denote the inverse Fourier-transform of y^2a;(/c)a(fe, a) by 'ijj{x, a). 
By the definition in (j^, we in fact then have il^{x,a) = J2y(p{x — y)qy + iapx- The 72- 
clustering assumption of the {q,p) variables then implies that the n:th cumulant l(yi = 
Q)k[^{x + yijiTi),... ,'ip{x + yn,o'n)] is square summable in y for every x. Thus one can 
safely take Fourier-transform with respect to the variables y, and the result is a function in 
L'^iry for every fixed x G Z. Since the second assumption implies that the function is slowly 
varying in x, on the scale we can thus conclude that for each a G {±1}"' there is a 
function Fn{x, k) such that it is slowly varying in the lattice position x on the scale it is 
L^-integrable in /c G T”, and 

K[y{ki,ai ),..., y{kn, an)] = Y ^^Fn{x, ki, k 2 ,..., kn, a). (40) 
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In fact, there are many such functions; by using the periodic lattice A as a middle step, we 
find that, for instance, any convex combination of the n functions, defined for Iq G {1,..., n} 
by 



j/SA'i 


(41) 


will work. 


In the translation invariant case, when the scale of spatial variation e —)> 0, is also 


independent of x and hence can be thought of as a proper function in L^(T"'). Since ^(/c, a) = 


^J2uj{k)a{k^ a), this yields the following structure for cumulants of the a-fields for translation 
invariant, ^ 2 -clustering, initial data: 



(42) 


In case ui has zeroes, we recall our regularized definition of the field a and conclude that 


Fn must then also be zero at such points. Let Wn{k,cr) denote the values from the product 
Fn{k,a) Y\£=i{‘2uj{ki))~^/‘^ . Thus, if uj{ki) = 0 for some then Wn{k,a) = 0 as well. 


Therefore, although the cumulants are singular, their singularity structure is simple, en¬ 
tirely encoded in the d-multiplier. In contrast, by the moments-to-cumulants formula, the 
moments satisfy for any sequence J of labels 



(43) 


This has ever more complicated singularity structure as the order of the moment is increased. 


(The above discussion can be made mathematically rigorous by replacing the infinite lattice 


by a periodic d-dimensional lattice. See [12] for details.) 

1.5 Green—Kubo formula 

The Green-Kubo formula gives an expression for the thermal conductivity if the system has 
normal heat conduction, i.e., diffusive energy transport. However, it can also be used to 
inspect anomalous transport since it always allows studying certain relaxation characteristics 
of perturbations of thermal states. 

Let us recall here the argument in the diffusive case, assuming that energy is the only 
relevant ergodic invariant. The basic assumption is that after a relatively short thermalization 
period the macroscopic energy density e(x, t) evolves according to the Fourier’s law 


dte = dx {D{e)dxe) . 


(44) 


This is a nonlinear diffusion equation with a diffusion “constant” D which depends on the 


value of energy density e. The steady states are then given by microcanonical ensembles 
whose only variable is the uniform energy density e. 


Consider next an infinite system with a steady state whose energy density is e. Perturb 


the state, locally near origin, so that the perturbation has finite energy. Then, for normally 
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conducting systems, one would expect that the energy density at “infinity” always remains 
equal to e and that the final state of an infinite system should be the same as the initial 
steady state, since for finite systems of volume V the final steady state should be of uniform 
energy density e + 0{y~^). This implies that there should be a time to after which the energy 
density is ever better approximated by a solution to the linear diffusion equation 

dte = D{e)d‘^e . (45) 


For a linear diffusion equation, it is possible to find out the diffusion constant from the leading 
evolution of the spatial spread of the perturbation. Namely, since then 


e{t + to,x) ~ / dy- 


1 


where f dx x‘^\e{to, x) — e\ < oo, it follows that 




(46) 


lim - dx x^{e{t, x) — e) = D{e)2SE , 
t^oo t I 


(47) 


where 5E denote the excess energy of the perturbation, which is conserved and hence 5E = 
f dx (e(t, x) — e) for all t > 0. 

Translated into the present lattice setting, one analogously arrives at the formula 


^lim - E x\[H^{t),Hom = D{e)2 , 


(48) 


where K[a, b] = {ab) — {a){b) = {{a— {a)){b— (b))) denotes the second cumulant, i.e., covariance, 
of the random variables a, b defined by solving the time-evolution using equilibrium initial 
data with mean energy e. Namely, if we divide the equation by {Hq), the sum on the right 
hand side corresponds to the excess energy from a perturbation of the initial measure yo to 
the probability measure -^^y/ro- The perturbation is then localized near the origin. Since the 
left hand side is a discrete version of the left hand side of (47), the limit in (48) should hold 


if the relaxation occurs via diffusion of energy, as given by (44). 


To get the standard Green-Kubo formula, we next also assume equivalence of the infinite 
volume microcanonical and canonical ensembles, and use here initial data distributed accord¬ 


ing to a canonical ensemble at a temperature T. Then e in (48) denotes the corresponding 


(uniform) equilibrium energy density, e = e{T) = Making this change of variables in 


(44) implies the standard Fourier’s law. Explicitly, we define the temperature distribution 


corresponding to a given e{x,t) by inverting the above function to yield T{t,x) = T{e{t,x)) 
which then satisfies 


de{T) 

dT 


dtT = {K{T)d^T) , 


(49) 


where the thermal conductivity is given by k{T) = D{e{T)) ^ . In the finite volume A 

canonical Gibbs state, we have e\{T) = {Ho)\ = Z~^ Jd^qd^pe~^^"^Ho where Z denotes 
the partition function, Z = f d^qd^pe~^^'^. Therefore, 


^ {{HHo)a - {H)a{Ho)a) = j;KA[i^.(0),i7o(0)] 




2^2 


(50) 
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Taking the infinite volume limit A —)> Z, and then using (48), thus yields 


i^iT) = ^ j x^K[H^{t), Ho{0)] . (51) 

xGZ 

The above argument shows that to find the numerical value of the thermal conductivity at 
a given temperature, it suffices to study the asymptotic increase of the spatial quadratic spread 
of the energy density, when starting the system with samples taken from the corresponding 
canonical Gibbs state. Explicitly, one should study the observable 

= (52) 


show that it is 0{t) for large t, and then solve the leading asymptotics. 

For this, it is convenient to rewrite S{t) in terms of current observables. Namely, for any 
J satisfying (32) we have for t > 0 and all x G Z 


Ha;{t) = Hx{0) + f ds {Jx-I,x{s) - Jx,x+i{s)) , (53) 

40 

Hx{-t) = Hx{d) “ y ds {Jx-I,x{s) - Jx,x+i{s)) . (54) 


Since the Gibbs state is invariant under both spatial translations and the Hamiltonian time- 
evolution, and these transform the energy observables as expected, we have here 


S{t) = Yx^^^\Ho{d),H_x{-t)] = Y^^^\Hx{-t),Hom = S{-t). (55) 

3iGZ fcGZ 


Hence, S{t) is symmetric, and 

2S{t) = 25(0) + fdsK[Jx-i,x{s) - Jx,x+i{s),HQ{d)] 

x & -^0 
/■O 

^ ^ ^ / ds K[Ja;_i^3;(s) (s), 1^0 (0)] 

x&l 

= 25(0) + J](x2 - (x - 1)2) fds k[Jx-iA^), {Ho{-s) - Ho{s))] , (56) 


xEZ 


where we have used the earlier proven simple rule for translation of the current observable. 
Therefore, by (53) and (54), we have 

5(t) = 5(0) - ^(2x - 1) f ds f drK[Jx-i,x{0),J-i,o{r)-Jo,i{r)] 

^ x& -^0 J-s 

= 5(0)-h / ds / dr ^k[Jj,_i,j;( 0), Jo,i(r)]. (57) 


10 J-s 


ieEZ 


We define the current-current correlation function, also called the Green-Kubo correlator, 
by the formula C{t,T) = ^['^ 2 :-i,a:(d))'^o,i(i)]- From the above properties, we obtain 


16 







c{-t) = = X]xez'^['^o,i(i);'^i-a:,2-x(0)] = C{t), and thus we find 

from 


s{t) = S{0) + 2 [ ds [ drC(r) = 5(0) + 2t [ dr (l- -) C(r). 
Jo Jo JO ^ ^ ^ 


(58) 


If /Q°°dr \ C{r]T)\ < oo, we can rely on dominated convergence in (58) and conclude that then, 
by (I^, 


k{T) = ^ dr C{r-,T). 


(59) 


This is the standard Green-Kubo formula, and thus, in great generality, the integrability 
of the correlation function C{t; T) is sufficient to imply normal thermal conduction at the 
temperature T. 

However, we can also conclude that, if fgdr (l — j) C{r;T) is not bounded in t, then the 
energy relaxation cannot be diffusive at temperature T. Hence, the Green-Kubo correlation 
function C{t]T) can also be used to prove superdiffusion of energy and to check at which 
time scale the energy spread occurs: if /gdr (l — C{r;T) = O(t^), with p > 0, then the 
spatial spread at time t is However, the spread from the observable S{t) can be 

somewhat misleading for systems which have polynomial decay or multiscale behaviour. In 
these cases, S{t) could be entirely dominated by tail behaviour, or just by one of the many 
scales, instead of describing the actual speed of spreading of local perturbations. We will find 


such an example from the kinetic theory of FPU-/? chains, discussed in Sec. 3.5 


Let us conclude the section by computing an expression for the current correlation C{t) 
for the particle chains. Some care needs to be taken with the sum over x and the principal 
value integral in the definition of the current observable. To control their behaviour, let us 
insert fjdk6{k — 'ffiki) into the integrals defining J in (33). Then for each t there is some 
function ft on the torus, such that 


«:[J3;_i,x(0), Jo,i(t)] = / dke 
Jt 

as a principal value integral around zero. Hence, 

1 


- 1 
i27rxfc 


- (l -hicot(7rfc)) ftik ), 


(60) 




- 1 ,*( 0 ), Jo,i(i)] = x/t(0) + i ^ 


^^^27rxk ftm fti k) 




(61) 


In fact, here /t(0) = 0 due to energy conservation, as discussed above. Also, for any function 
F we have dj.j^„dk5{k - Yl,e.k()F{k)\-^^^ = J^„dk 8{Yftke)dk^F{k). 

Therefore, for the particle chains the Green-Kubo correlation function is given by 

^K[Jj;_i,a;(0), Jo,i(t)] = ^ —'^Wnicr)wn'{cr')[ dfc 5 ( ^ ] [ dk' 

XGZ n',n=2 a'\7^1 ) 

n' 

X -^^n'{k')(-l -|-icot(7ry~]fc£( 


e=i 


X K 


n n 

dki ^>n(/c)5s(]^a5(A:£,c7£)^^_^ , cr^)) 


e=i 


e=i 


(62) 
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where the weight Wn{o') = 1, for a G {±1}*^, unless n = 2 and ui = a 2 , in which case it is 
zero. By translation invariance, the final cumulant is proportional to 

and thus the integral over k' has to be treated similarly to what has been done above for the 
integral over k. Following the same steps as above, we thus arrive at the simplihed formula 


Jo,.(()]= E E E Wn{(T)Wn'{cT') 

n',n=2 o-'g{±l}n' o-e{±l}" 


x£jj 


X / dkdi'^ki] / dk'K Jnik,a]0),Jn'ik' 


(63) 




where for k G T”, a G {±1}”, we have defined 


Jn{k,a;t) = 

n 2 tt 


n 

^n{k)dt(j^at{ki, 

i=i 


<yt 


(64) 


Assuming that the harmonic terms dominate over the anharmonic ones, we can obtain 
the leading harmonic part of the above expression by considering only the term n' = 2 = n 
above and approximating dtat{k,a) « —iau}{k)at{k,a). This yields an approximation 


J2{k,a;t) Ri ^n{k){criu{ki) + cr2u;{k2))at{ki,ai)at{k2,a2] 

47r L 


(65) 


Since these terms appear in the integrand only when a 2 = —(7i and k 2 = —ki, for which 
cJia;(A:i) + a 2 Uj{k 2 ) = 0, only the term where the factor criu;(A:i) + a 2 Uj{k 2 ) is differentiated 
yields a non-zero contribution to the integral. Since then also ^ 2 {k) = oj{ki), it is possible to 
substitute the above harmonic current term by 


w(A;i) 


dvr 


U}'{ki)at{ki,ai)at{k2,cr2) 


( 66 ) 


in the integrand. Thus the harmonic part of the Green-Kubo correlation function evaluates 
to 


cr a 


E 4 

cr',(Tg{±l} 


^ / dk 


f dk' ‘^^-^uj{k) ^ uj{k'i)K ao{k, a)ao{-k, -a), at{k[,a')at{k 2 , -a') 

/'j’2 ZTT ZTT 


(67) 


As ijj'{—k) = —uj'{k), we can simplify the expression (67) further to 

,,J{k) 


dk / dA:' 
/t 2 


/t 


27r 


-uj{k) 


27r 


io{k' 2 )n aoi-k,-l)ao{k, l),atik[,-l)at{k 2 ,1 ) 


( 68 ) 


Here u:{k) = {k) which is equal to sin(27rA:) for chains with stable nearest neigh¬ 

bour interactions, using the parametrization given in Section 1.2 Thus the harmonic in¬ 
teraction part of the Green-Kubo correlation function for nearest neighbour chains is equal 
to 


—4 r2 
iO 0 K 


L4t 


dA:sin(27rA;)ao(—A:, —l)ao(A:, 1), / dfc'sin(27rA;2)at(A;(^, —l)at(A; 2 ,1) 


/T 2 


(69) 
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This result coincides with the formula used in Section 3 in |T6]. As explained there, the 
above time-correlator can also be computed by solving the evolution of the Wigner function 
starting from a suitably perturbed initial Gibbs state. Namely, consider the Gibbs states with 
perturbed Hamiltonians H^{q,p) = H{q,p) — TeJ[q,p) where 


= / dA:sin(27r/c)oo(—A:, —l)ao(/c, 1) (70) 

Jt 

and denote the expectation over the canonical Gibbs state at temperature T and the Hamil¬ 
tonian by (•)^^^. For instance by using approximations by finite periodic lattices, it follows 
that 


de { I dk'sm{27rk2)at{k'i,—l)at{k2,l) 


(e) 


£ = 0 




' dA:sin(27r/c)ao(—A:, —l)ao(A:, 1), / dA:'sin(27rA:2)at(A:(, —l)at(A: 2 ,1) 
L./T Jt^ 


(71) 


Note that the time-evolution in both expectations is the same, determined by the original 
Hamiltonian H. Hence, computing the left hand side involves solving the original problem 
with non-stationary initial data. The perturbation in the initial data can be captured by 
changing /3uj{k) to I3uj{k) — sa sm{27rk) in the harmonic part multiplying ao{k, —a)ao{k, a). 
In particular, also the perturbed initial state is translation invariant in space. 


2 Phonon Boltzmann equation of spatially homogeneous an- 
harmonic chains 

2.1 Identifying the collision operator 

The goal of the section is to derive the form of the phonon Boltzmann collision operators for 
the particle chains. We do this only for spatially homogeneous initial data which suffices to 
cover the evolution of the Green-Kubo correlation function, as explained in Sec. 1.5 Allowing 


for inhomogeneous initial data would be an important extension, but the inhomogeneity 
introduces new multiscale effects which we wish to avoid in the present contribution. See 
the discussion in Section 6 of m for an example of the problems associated with spatially 
inhomogeneous initial data for a nonlinear evolution equation similar to the FPU chains. 

Let us first recall the parametrization of the spatially homogeneous correlation functions 


Wnik,a-,t) satisfying (42), 


n 

K[atiki,(Ti),... ,at{kn,(Tn)] = 5(^'^ki^Wnik,a]t). 


(72) 


e=i 


As discussed earlier, this requirement does not fix the functions Wn completely, and here we 
resolve the ambiguity by requiring that they are constant in the first variable, ki. This is 
equivalent to using the definition 


Wn{k,a;t)= / dk[K[at{k[,ai),atik2,a2),...,at{kn,(rn)]- 


(73) 
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We can then also obtain the Wigner function W{k]t) of the state, as used in earlier works, 
by setting W{k; t) = W 2 {{k', k), (—1,1); t) for all /c', /c G T. 

The time-evolution of the correlation functions Wn can now be solved from (37). Recalling 
the earlier notations involving a sequence J of labels, we have 


dtK[at{k, a)j] = ag) :at{k, . 

i£j 


Since the derivative of the phonon field satisfies (23), we thus obtain 


dtK[at{k,a)j] = -iO„(A:, u) jK[at(/c, cr) j] 

Z [ dk'6ike-k[-k'^)<^3{ke,k') 


— 1 


e&J c.'e{±i}2' 

X {at{k[, a[)at{k 2 , ( 12 ) :at{k,ay'^^:) 


— 1 


E L dk'6ike-^y)^yki,k^) 


e&J (T'e{±i}3 i=i 

X {at{k[,a[)at{k2^cr2)at{k'^,ay :at{k,ay'^^:) , 


(74) 


(75) 


where n = \J\ and we have used (36) to simplify the linear, free evolution, term into a 
multiplication by the function 


0.n{k,a) = y^^aiu{ke) 


(76) 


e=i 


The remaining expectation values can be represented in terms of cumulants using (35). To 
derive the Boltzmann equation, we will need only the first four cumulants but let us postpone 
the detailed expansion for later. 

We next exponentiate the free evolution part by integrating the identity 


= e-^^rr{k,a)j{t-s) a)jK.[as{k, a)j] + dsK[as{k, u) j]) . 


This results in the following Duhamel perturbation formula 

rt 


iy^ae [ ds e-iW(fc,<x)Ri-.) ^ f _ A;')$3(fe^,fc') 

/P.7 ^'c7-l-n2 


ieJ <T'e{±i}2 

X {asik[, a[)asik2, :as{k, ay^^:) 

i ^ ap fds e-iW(fc,<x)Ri-.) ^ f ^ A:' 

teJ (T'e{±i}3 i=i 

X {as{k[, a[)as{k2, a'2)as{k'^,a'y :as{k, ay\^:). 


)Mke,k') 


(77) 


(78) 


The case with n = 2 and J = ((fci, —1), (^ 2 ,1)) determines the Wigner function W{k;t) 
after integration over kp. However, in this case Hn(fe, a)j = 07 (^ 2 ) — io{ki) which is identically 
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zero whenever ki + k 2 = 0. By translation invariance, all three terms on the right hand side 
of (78) are now proportional to 5{ki + k 2 ), the last two since the remaining expectations are 
proportional to 6 {ki/ + K) with [(.') = J\£. Therefore, all the phase factors involving Qn 
are equal to one in this case. 

Let us next assume that the fields have been centred, i.e., that {at{k, a)) = 0. In general, 
for a spatially homogeneous case we could have {at{k,a)) = 6 {k)ct{o'). However, for instance 
in the FPU model we have normalized the fields so that at{0,a) = 0 for all t, which implies 
Ci(cj) = 0. Under this simplifying assumption, we obtain using ( |35[ ) and the obvious shorthand 
notations for the fields that 

{a[a 2 :a:) = 4, a] , (79) 

(a'^02®3 -O--) = ®] ) (80) 


where a = as{ki/,a£/). Therefore, integration of (78) over ki and writing k = k 2 yields 
W{k;t) = W{k;0) 

-i [ ds V 0-3 / dk'S{k'^-a 3 k)dC^k'i)^ 3 {k')W 3 {k',a'-,s) 

Jo 7t3 ^ 


a'e{±i}3 


2=1 

4 


— i ds 1 T 4 / dk' 6 {k'^ —a 4 ^k) 6 (S^ k'i)^ 4 {k')W 4 {k',a';s) 

Jo Jt^ 


'e{±i}4 


- 3i / ds 


E 

'e{±i}4 


i=l 

4 


Ij4 


dk' 5{k'^ - (j'^k)5C^ K)^{k'i + A:^)^4(fe') 


2=1 


X IU2((fe'l, k' 2 ), (o-'i, 0-2); s)W 2 {{k'^, k' 4 ), (0-3, C74); s ), 


(81) 


where we have enumerated the sum over i £ J by using a'j^ = —ai, when clearly ki = —cr'^k 
and ki> = a'^k. We have also used here the permutation invariance of and the obvious 
antisymmetries of d>„, n = 3,4, under a reversal of the sign of the first component. 

The last term in (81) can be simplified further using symmetries of <^4 to 

-3i 


ds cr 2 lT 2 ((-cr 2 A:,cr 2 /c), ((Ti,o- 2 ); s) 

(Te{±i}2 


X ljlk'M-k',k',-k,k) W 2 {i-k',k'),{a[,a' 2 y,s). 

<T'e{±i}2 


It 


(82) 


Inserting the definition of the FPU interaction amplitude <^4 given in (27) then yields 

rt 


i 6 A 4 a; ^ / ds | sin(7r/c) | E 0'2W2((-fT2fe,CT2fe), (cri,cr2);s) 

(Te{±i }2 


10 


/d/c'I sin(7rfc')| V W 2 {{-k', k'), {a[, a'^)] s). 
JT 


(83) 


' 6 {± 1 } 


Since now uj\ sin( 7 rA;)| = oj{k)ly /2 and 


a 2 W 2 {{-<T 2 k,(J 2 k), (cri,cr 2 );s) = W 2 {{-k,k), (1, l);s) - W 2 {{-k,k), (-1,-1); s), 

7e{±i}2 
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for the FPU chains we may write this term as 


i3A4a; ^ / dsLj{k) aW 2 {{—k,k),{a,a)]s) 

X fdk'uik') V W2ii-k',k'),{a[,a'^y,s). 
CT'e{±i }2 


( 84 ) 


It is customary in kinetic theory to ignore the term (84) in the Boltzmann equation. At 


the first glance, it is not obvious why this should be so: the term does not involve any non- 
Gaussian factors and it has an apparent magnitude 0{X/^t). Although no full mathematical 
study has been made about the effect of the term on the kinetic time scales 0 (AJ ), the 
following argument suggests that, indeed, it behaves better than suggested by its apparent 
magnitude which is O(AJ^) on the kinetic time scale. In fact, the argument indicates that 
the term is 0 (A 4 ) uniformly in t and k, and hence could be safely neglected when A 4 <C 1 . 


We study the two kF 2 -dependent factors in (84) separately. The second factor is 


dkuj{k) ^ W2{{-k,k),{ai,a2);s) 

^ CTe{±i }2 

= / d'^k y/u!{ki)uj{k2) ^ K[as{ki,ai),as{k2,(J2)] 


/T 2 


o-e{±i}2 


= 2 / d'^kuj{ki)u){k2){qis,ki)q{s,k2)), 


(85) 


where we have applied the inversion formula given in ([^ and the assumption that {as{k, a)) = 
0. This is equal to 2{{uj * g(s))(0)^) which for a finite periodic lattice would correspond 
to 2|A|“^X)x((^ * 9('S))(^)^) = 2|A|“^(X)2,_y g'a:(s)a(x - y)qy{s)), by translation invariance. 
Therefore, the second factor is proportional to the average harmonic potential energy, and 
we would expect it to be slowly varying and to equilibrate to the corresponding equilibrium 
value as s —>■ 00 . In particular, it can be bounded by the total energy density which is a 
constant, so the value of this term is 0 ( 1 ) uniformly in s. 

The first factor in (84) depends only on the “off-diagonal” elements of the second order 
correlation matrix. By (78), they have rapidly oscillating phase factors e“'^^ where U = 


2auj{k). If we assume that this oscillatory behaviour continues for all s, it seems reasonable 


to expect that the value of (84) can be well approximated by using this oscillatory factor and 


assuming that the remainder is so slowly varying that it can be replaced by a constant. This 
yields a factor which is 0 (A 4 ) times the integral U Jg*dse“'^^ = i(l — e“'*^) which is bounded 
by 2 for all t and k. Of course, the replacement of a slowly varying term by a constant is not 
totally accurate here, but it seems reasonable to assume that under the present conditions 
the resulting corrections would also be 0 (A 4 ), at least up to the kinetic time scales. 

The rest of the argument leading to the phonon Boltzmann collision operator is fairly 
standard. The clusterings relevant to the third cumulant are 


{a'ia2 :aia2'.) = K[a[,a2,ai,a2] + K[a'i,ai]K[a2,a2] -|-a2]K[a2, ai] 


( 86 ) 


and ( 0(^0203 :aia 2 :) whose expansion has only terms containing higher order cumulant factors. 
We continue to assume that the oscillations of the higher order cumulants are sufficient to 
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suppress any term which involves integrals of Wn, n > 2, over any of its arguments. The 
clusterings relevant to the fourth cumulant are (a'^a2 :aia2a3:), which produces only terms 
containing higher order cumulant factors, and 

(a'l^ag :aia2a3:) = ^ «i]^K(2)> «2]^tK(3), as] + (higher order terms), (87) 

7rSPerm(3) 


where the sum goes through all permutations vr of three elements. 

In summary, we then obtain the following approximations from the third and fourth 
cumulants. 




-2i^a£ dse ^)^)3(fc)]J ^ W 2 {{-ki,ki),{a ,ai)-,s), (88) 

1=1 ijiea'&{±l} 

W4{k, (j; t) ^ u; 0) 

,cji);s), (89) 

t7'e{±i} 


4 .t 

+ 6iJ]^7W dse-'''4('^''^)(*-^)<I>4(A:)n ^ IT2((-fei, (a', < 

^=1 i^^a'&l±l} 


where we have simplified the result by applying the invariance properties of the amplitude 
factors d>3 and <^4. We insert these approximations in (81). The terms depending on the 
cumulants at time zero are again assumed to be negligible since they involve highly oscillatory 
integrands. 

Since ^2 ((—/ci), (u', Uj); s) is equal to W{aiki,s) plus a rapidly oscillatory 

(cTj, fTi)-term, this results in the approximation 


W{k-t) ^ W{k-,Q) 

+ 2 f ds Y, ^ 3 / dA:'(5(A:^-U 3 A:)( 5 (^A:')(-^' 3 (A:')^) [ dt'e-‘^3(fc'.a)(i'-s) 


(7e{±i}3 


i=l 


2 = 1 


+ 6 / ds V (74 / dk' 5{k'^- aAk)5{Y k'i)^A{k'? [ 

Jo S Js 


7e{±i} 

4 

X ^CJ£]JlT((Ti/c',s) , 
e=i i^e 


(90) 


where we have used Fubini’s theorem to exchange the order of the time-integrals. If we swap 
the sign of both a and k', clearly also Qn{k', a), n = 3,4, change their sign. Therefore, in the 
above formula, we may replace both of the terms g-i(('-s)Q 1 ~ 

for large t. 
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Inserting the explicit expressions for thus results in the formula 
W{k]t) « VF(fc;0) 


2r^-3 


“t“ 2 2 TTAgCU 


(7e{±l}3 i=l i=l 


aiU}{k'i)) 


]J I sin(7r/c^) | ^ ctsu^ ]J Vl^(A:',, 


£=i 


£=i 


i-^l 


t 4 4 

+ 24:7r\‘llj~^ f ds ^ f dk' 5{k'^ —k)5(^^a[k[)5('^^aiUj{k[)) 

fTe{±i}4 j=i i=i 

4 4 

X ]J|sin(7r/c^)|^cj4fT£]Jiy(A:',s). (91) 

£=i £=i 

The right hand side is an integrated version of the solution to the homogeneous Boltzmann 
equation 

dtW{k-t) = AiC3[IT(-;t)](fc) + AlC4[IT(-;t)](fc), (92) 

where, using shorthand notations uji = uj{k£), IT^ = IT (ki), ^ = 0,1, 2, 3,4, 

Cs[W]{ko) = / dk 6{ko + aiki + (72k2)S{u;o + criui + a2Ui2) 

_.7r:i,4T2 


<7e{±i}2' 

2 

X I sin(7rA:£)| x (IT 1 IT 2 + cJiIToVh '2 + trsIToITi) , 


(93) 


e=o 


and 


„ 3 3 3 

C4[IT](A;o) = 487ra;“'^ ^ / dk6{ko + 'y^^(Tiki)5{uJo + y^Q-jCJi) J] \ sin(7rfc£)| 

f7e{±i}3*=i *=i f=o 

X (IT 1 IT 2 IT 3 + CTiIToW'2lT3 + CJ 2 IT 0 IT 1 IT 3 + a3lToITiIT2) . (94) 

Other interactions can be treated similarly. For instance, the onsite perturbations have 
a: 3 {ki,k 2 ) = A 3 and a/i{ki, k 2 , k^) = A 4 and, for general dispersion relations to, yield the 
collision operators 

r 

C3[f^](^o) = 47r ^ / dk5{ko + aiki+a2k2)d{uo + aiU}i+a2UJ2)Y\_2 — 

tTe{±i }2 e=o 

X (IT 1 IT 2 + aiIToIT2 + fT2VFoITi) , (95) 

and 

3 3 3 

C4[IT](A:o) = 127r ^ f dkSjkp + y^^criki)6{uJo + y^u^cui) J] ^ 

(Te{±l}3 i=l i=l £=o 

X (ITiIT2hF3 + fTiIToIT2lT3 + (T 2 VF 0 IT 1 IT 3 + (T 3 VF 0 IT 1 IT 2 ) . (96) 

Detailed derivations of these onsite collision operators using standard perturbation expansions 
can also be found in the literature: the third order collision operator in [T] and the fourth 
order collision operator in [16] , 
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2.2 Solution of the collisional constraints 

The above derivation of the Boltzmann collision operator is not mathematically rigorous. 
In fact, the argument used for neglecting the higher order terms and replacing “t — s” by 
“oo” in the derivation of the energy conservation (5-function are at present mathematically 
uncontrollable approximations. Here we assume that we are working in a regime in which these 
terms indeed can be neglected. Hence, instead of trying to find out under which conditions 
this would work for our particle chains, we start from the homogeneous Boltzmann equation 
and proceed to study its solutions to obtain predictions which can later be compared with 
other results on the chains, such as computer simulations. 

The weakest part of the argument lies in the rule that “any integral over one of the 
fc-components of the rapidly oscillating factor leads to integrable decay in t and 

hence can be neglected”. This is particularly suspect if the integration is performed over a 
one-dimensional space; typically, even the best decay estimates over a (i-dimensional torus 
yield Jrjp^dA: which is not integrable ii d = 1. The decay of such os¬ 
cillatory integrals is closely linked to the behaviour of the frequency function around 

its zero set; consider, for example, the explicitly integrable example JqcIs f = 

I Q — 1) which is uniformly bounded in t, if is integrable, but which 

can blow up as fast as 0{t), for instance, if the set of k with il.n{k) = 0 has a non-zero 
measure. 

These subtleties are also reflected in the Boltzmann collision operator: the two (5-functions 
need to be carefully integrated over before one can use the equation. In the well-studied 
continuum setup of Boltzmann equations, the dispersion relation is given by k G 
and the integral contains S{ko + ki — k 2 — /c3)(5((/cg -|- /cf — — kf)/2) which can be explicitly 

integrated over, yielding the “standard” rarefied gas collision operator. 

For a lattice system, the solutions to ^n{k) = 0 are not so easy to handle. As we will 
show below, there are however explicit solutions for the one-dimensional nearest neighbour 
dispersion relation, relevant to the FPU chains. After the solution manifold has been found, 
one still needs to choose some parametrization of the manifold to integrate over the energy 
conservation (5-function, since this will yield an additional factor to the remaining integrand. 

In fact, even at this part of the procedure, some care is needed since the trivial solutions 
to the collisional constraints, to be discussed later, do not contribute to the total collision 
operator but would result in infinite factors if integrated over using the above rule. Hence, 
the trivial solutions need to be removed from the solution manifold before integrating out the 
collisional constraint (5-functions. 

We begin with a result which shows that nearest neighbour dispersion relations are, in 
fact, somewhat pathological: they suppress all collisions involving three phonons. 

2.2.1 Cs = 0 for nearest neighbour dispersion relations 

Consider the nearest neighbour dispersion relation, Lo{k) = (1 — 2(5cos(27rA:))^/^ with 0 < (5 < | 
(we ignore the overall prefactor here, since it does not affect the solution manifold). To esti¬ 
mate the collision energy U3 = ujQ+aiUJi+a 20 J 2 -, consider the following parametrization of oj as 
the magnitude of a complex number: Since 0 < 2(5 < 1, we can define r = 2ai'cosh(2(5)“^ > 0. 
Then by explicit computation 

uj{k) = \z{k)\^ for 2:(/c; (5) = \/5(e’’— . (97) 
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The above parametrization and the triangle inequality yield an upper bound for differences 
of dispersion relations. Namely, for k,q G T, 

\oj{k) — uj{k + g)| = I \z{k)\ — \z{k + (?)| | < \z{k) — z{k + (?)!= e“^\/5|l — . (98) 

Since |1 — = 2 — 2cos(27rg) = 5~^{25 — 25cos(27rg)) < 5~^u:{qY, then 

\uj{k) — uj{k + q)\ <Q~'^oj[q). (99) 

If 5 < ^, we have r > 0, and hence e“^ < 1. In addition, then Ljj{k) > c^, with = 
(1 — 25)^/^ > 0. Thus if the momentum constraint k^ + aiki + (72^2 = 0 holds, we have 


Ha 

= ojQ + uji + 0 J 2 > 3c5 , 

for 

fji = 1 = fT 2 , 

( 100 ) 

ria 

= cjo - - 1 - a ;2 > (1 - e“^)c 5 , 

for 

- fji = 1 = 0-2 , 

( 101 ) 

ria 

= Wo - 1 - - 60)2 > (1 - e“'’)c 5 , 

for 

(7l = 1 = -CJ 2 , 

( 102 ) 

ria 

= Wo - wi - W 2 < -(1 - e“'’)c 5 , 

for 

(7l = -1 = (72 . 

(103) 


Therefore, with a nearest neighbour interactions and pinning, the momentum constraint keeps 
the three phonon energy well separated from zero, at least a distance (1 — e~'")cs apart. In 
such a case, one should set C3 = 0 in the phonon Boltzmann equation. 

If 5 = 2, such as in the FPU models, we have r = 0 and cs = 0, and equation (99) does 
not rule out the existence of solutions to both constraints. However, it is possible to improve 
the bound (98), by recalling that an equality in the triangle inequality holds if and only if one 
of the numbers is a non-negative multiple of the other one. This requires that for the given 
k, q there should be i? > 0 such that (1 — = R{1 — , where ki = k, k 2 = k + q, 

or vice versa. For instance by a geometric argument, it is straightforward to see that this 
happens if and only if /ci = 0 or k 2 = 0. 

As above, by inspecting the four sign combinations, we thus find that IIIsI > 0 unless 
/co = 0, fei = 0, or k 2 = 0. If fco = 0, we have fis = 0 identically for cti = — 1T2, and only at 
the point fci = 0 = A;2, if ui = cr2. If /cq / 0, we have Us > 0 for ai = 1 = <72, and fla = 0, 
only if one of the following conditions is satisfied: 


1. —(7i = 1 = (72, and k 2 = 0, 

2. (7i = I = —(72, and ki = 0, 

3. (7i = — 1 = (72, and /si = 0 or ^2 = 0. 


However, these solutions do not contribute to the collision operator in the FPU models as 
the collision operator C3 defined in (93) contains a prefactor nr=o I sin(7rA;£)| which is zero on 
any of the above solution manifolds. Therefore, for the FPU-chains, the lowest order kinetic 
theory would imply using C3 = 0. 

The above argument straightforwardly generalizes to higher dimensions. As shown in 
Appendix 18.1 of [T], also higher dimensional nearest neighbour dispersion relations with 
pinning have no solutions to the constraints in the three-phonon collision operator and, hence, 
C3 = 0. However, this property depends on the dispersion relation. For instance, consider 
the next-to-nearest neighbour interaction with a{k) = (1 — cos(27rA;))^. Then the dispersion 
relation is uj{k) = 1 — cos(27rA:) and for ai = 1 = —(72 we thus have 


Q 3 = ujQ LOi — UJ 2 = I — cos(27r/co) — cos(27rA:i) -|- cos(27r(A:o -|- ki )), (104) 

which is zero ki = ^ ~ ^0- This introduces solutions to the three-phonon collisional con¬ 
straints in the next-to-nearest neighbour case. 
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2.2.2 C4 for nearest neighbour dispersion relations 

From now on, we only consider the nearest neighbour dispersion relations. We do not claim 
that this case is representative of the general case: as the previous examples indicate, other 
types of behaviour might appear for other dispersion relations. 

Consider thus the collisional constraints in the four-phonon collision operator, assuming 
a nearest neighbour dispersion relation. If iTj = 1 for all i, by the above estimates, we have 
^4 > 0, unless 5 = 5 and ki = 0 for every i = 0,1, 2, 3. 

If Uj = —1 for all i, we resolve the momentum constraint by integrating out which 

yields k^ = ko — ki — k 2 - Then 

Q4 = LOq — iOl — U!2 — CO3 

= ^Isiko, ki, ko- ki] 1, -1, -1) Qsiko - ki, ^ 2 , ko - ki - k2] 1, -1, -1) . (105) 

By the estimates derived in section |2.2.1t we then have 114 < 0 uniformly in the case with 

pinning. In addition, if 5 = ^, the only way to get 114 = 0 is that both of the terms above 
are zero, which happens only if fei = 0 or ki = ko, and k 2 = 0 01 k 2 = ko — ki, i.e., k^ = 0. 
This implies that the only solutions to 04 = 0 are the trivial solutions where two of the wave 
numbers ki,k 2 , k^ are zero and the remaining one is equal to ko. 

The final degenerate case is found if = 1- Then only one of cTj is negative; 

for notational simplicity let us suppose it is cii (the other cases can then be obtained by 
permutation of the indices). The momentum conservation implies then that ki — ko = /c 2 + ^ 3 ; 
and hence 


124 — ujo — T ^2 T W3 

= 113 (^ 0 , ki, ko-ki] 1, -1,1) -h 03 (^ 2 , k3, k2 + ko] 1,1, -1) 


(106) 


Thus ^4 = 0 only for the trivial solution in the unpinned case, namely, if ki = ko and 
A:2 = 0 = ko. 

The remaining sign-combinations have 1 + X]i=i = 0) i-®-) l^e polarization of phonons is 
preserved in the collision. These are also called collisions which conserve the phonon number. 
From the above results we can already conclude that every term in the collision operator 
which is not phonon number conserving has only trivial solutions to the energy constraint, 
and these solutions all require that two of the ki, i = 1,2,3, are equal to zero and the last 
one is equal to /co¬ 
in particular, there are no trivial solutions in the pinned case, so these can again be 
safely neglected. The same holds for the FPU chains, due to the prefactor n£=o I sin(7r/c£)|, 
and since the solution involves a two-dimensional integral and an integrand which has only a 
point-singularity. To have more convincing an argument, one should choose a regularization 
of the (5-function and to show that the contribution from the resulting ordinary integrals 
vanishes as the regularization is removed. An example of such a procedure is given in m 
Section 3] where it is used for showing that the trivial solutions in the number conserving 
case do not contribute to the linearized collision operator. 

Here we are interested in the collision operators defined in (94) and (96). Using the 


permutation properties of the integrands multiplying the constraint (5-functions, we thus arrive 
at the following simplified forms for the collision operators: in the FPU chains, the standard 
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kinetic argument yields 


C[kk](/co) = Ott (2A4)^w ® / dk 5 {kQ + ki — k2 — k^) 5 {ujQ + oji — UJ2 — 0 Jz)W oji 

Jt3 


e=o 


X {W1W2W3 + W0W2W3 - W0W1W3 - W0W1W2) , 
and for onsite nonlinear interactions 


( 107 ) 


Ovr 


C[W]{ko) = / d/c(5(/i:o + ki - k 2 - k3)S{oJo + uii - UJ 2 - 003 ) ^ 


4 Jj3 

X {W1W2W3 + W0W2W3 - Ik^oVkiVks - W0W1W2) . 


£=0 


(108) 


The two operators are nearly identical, differing only by the powers of the uj{ki) factors. 
This difference however has an important consequence to the properties of the solutions to 


these two equations; the onsite equation (108) predicts finite thermal conductivity from the 
Green-Kubo formula, whereas the FPU chain equation predicts an infinite result. Qualita¬ 
tively, one can understand this by noticing that in the unpinned case w has zeroes which leads 
to enhancement of collisions for the onsite anharmonicity, but which suppresses them in the 
FPU case. However, the qualitative argument is not sufficient to determine the magnitude of 
the effect, and we need to do a more careful study of the properties of the solutions to obtain 
predictions about the thermal conductivities. 


2 . 2.3 Integration of the collisional constraints in FPU chains 

Let us begin with the unpinned nearest neighbour models relevant to the FPU chains. Since 
then u!{k) = \/2w| sin(7rA:)|, it is analytically simpler to reparametrize the wave-number inte¬ 
grals by changing variables from k € T to p = 27 rk which belongs to the half-open interval 
I = [0,27r). Then uj{k) is replaced by tD(p) = \/2sin(p/2) where the absolute value is not 
needed since it was possible to restrict the values of p to an interval where sin(p/2) is positive. 

Although very useful for derivation of explicit parametrizations of the solution manifold, 
it should be stressed that some care is needed with this choice: it is crucial below that all 
arithmetic involving p € I is performed modulo I. For instance, if pi = tt/2 and p 2 = 37r/2, 
we have to use pi — P 2 = +'k in uj{pi — P 2 ) to get its value correctly. This rule of “modulo I 
arithmetic” will be applied without further mention below. 

We make the change of variables also in the spatially homogeneous Wigner function and 
consider Wt{p) = lUt(p/(27r)). Let us for simplicity drop the tilde from here and from the 
reparametrized dispersion relation, and denote them simply by Wt{p) and uj{p). With these 
conventions, the phonon Boltzmann equation of the FPU chain becomes dtWt{p) = C[Wt]{p) 
where 

3 

C[W]{po) = [ dpidp2 6{n{p))T\uji 

^ 1=0 

X {W1W2W3 + W0W2W3 - W0W1W3 - W0W1W2) , (109) 

with Ldi = uj{p£), Wi = W{p£), and 

Q{p) = ujo + uji - UJ2 - u}3 , p3 = {po + Pi - P2) mod I. (110) 
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The detailed solution of the collisional constraints, n(p) = 0 in (110), can be found in 
m- By Corollary 3.3. there, we can conclude that both constraints are satisfied exactly for 
those Pi, = 0,1, 2, 3, which satisfy one of the following three relations 


1- P 2 =Po, 

2- Pi = P 2 , or 

3. Pi = h{po,p 2 ) mod I, where h is defined using the standard, non-periodic, arithmetic in 

,, , v — x ( \y — x\ y + x\ 

h{x,y) = — - -h 2 arcsm( tan J^ cos —j (m) 


with arcsin denoting the principal branch with values in [—7r/2,7r/2]. 

The solutions satisfying item 1 or 2 are called perturbative or trivial, while the solution 
satisfying 3 is called non-perturbative. The nomenclature can be understood by expanding 
the constraint Q around small values of G M using (-o{p) ~ 2“2 |p| (1 — p"^/24:). The resulting 
equation then has only 1 and 2 as its solutions. 

The remaining energy conservation (5-function can then be formally resolved by integrating 
over a suitably chosen direction in the pi,p 2 -variables. For instance, choosing the pi-integral 
for this purpose would yield for any p 2 / po and for any continuous periodic function G, 


dpi6{n{p))G{po,Pi,P2) 


1 


\d2^iP0,P2,P2)\ 


G{P 0 ,P 2 ,P 2 ) + 


1 


\d2^{po,h{po,P2),P2)\ 


G{po,h{po,p 2 ),P 2 )- 


( 112 ) 


However, this procedure needs to be used with some care: for instance, in the present case 


fjdp 2 \d 2 ^{po,P 2 ,P 2 )\~^ = oo, and thus the first term on the right hand side of (112) would 
typically diverge when integrated over p 2 . 

This problem is resolved in the collision operator by the alternating signs in its integrand 
which guarantee that the integrand vanishes at the trivial solutions. This can even be proven 
rigorously if 1/W is sufficiently regular, say twice continuously differentiable, by using the 
following equivalent form for the collision operator: 


9 


TT 


/ dpidp2 5{Q{p)) \{{uiWi) (lFo“' + - ^2~' - ■ (113) 

e=o 


However, we do not wish to go into more detail here, but just add the contribution from the 
trivial solutions to the class of terms which are assumed to be negligible in the kinetic theory 
of FPU lattices. 

If the trivial solutions are neglected, we can use results from HZl, which rely on the explicit 
form of the non-trivial solution h, and obtain the fully integrated form 


C[W]{po) 



1 

V^F+(po,P2) 


3 


e=o 


X {W 1 W 2 W 3 + W 0 W 2 W 3 - IT 0 VF 1 IU 3 - W 0 W 1 W 2 ) , 


(114) 
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where pi = h{po,p 2 ), P 3 = Po + Pi - P 2 , and 


^ / X y\2 ... X . 

x± (ic, y) = ^^cos — + cos -j ±4 sin — sin 


(115) 


(See for instance, Lemma 3.4 in |17) for more details. From the Lemma a weak convergence 
of the integrals with a regularized J-function <5e(n) = e 7 r“^(e^ + e > 0 , and assuming 

continuity of W, can be established.) 

The function FL dehned in (115) is related to the change of variables which corresponds 
to using p 2 instead of pi to integrate out the energy 5-function. Namely, as proven in Lemma 
3.5 in dZ], 


/ dpodp2^=^== 

112 y'F+{po,P2) 


Gipo,h{po,P 2 )) = 2 [ dpodpi ^^^^^^==^=^G(po,Fi) , (116) 

Jf y'F_{po,pi) 


for any G for which either of the two integrals converges. The characteristic function restricts 
the integral to the subset of in which the argument of the square root is positive. Let 
us also mention that the change of variables becomes more involved if G also depends on p 2 
directly since each pair {po,pi), for which F_(po,Fi) > 0 , has two distinct values p 2 solving 
the energy constraint. Further details about these solutions can be found from the proof of 
the above mentioned Lemma 3.5. 


2.2.4 Integration of the collisional constraints for onsite nonlinearity 

To allow a comparison, let us also consider the case with an onsite nonlinearity. If the 
harmonic part has no pinning, then the above FPU discussion applies immediately, since 
these models have the same dispersion relation and differ only by the factors in the integrand. 
We thus find the following integrated form of the collision operator, if 5 = ^ in the general 
nearest neighbour dispersion relation: 


C[W]{po) 



1 

\/^+(P0,P2) 


3 




X (IU 1 IU 2 IU 3 + W 0 W 2 W 3 - W 0 W 1 W 3 - VF 0 IU 1 IU 2 ) , 


(117) 


where pi = h{pQ,p 2 ), P 3 = Po + Pi — P 2 , and F+ has been defined in ( 115[ ). (We continue to 
neglect contributions from the trivial solutions.) 

The analytic structure of the solution manifold gets more complicated once the dispersion 
relation has pinning, i.e., when 5 < Nevertheless, it is still possible to find a function 
h{pQ,p 2 ) such that the non-perturbative solution is parametrized by a condition 


Pi = Hpo,P2;S) . 


(118) 


In other words, the enumeration by the three conditions in the previous section continues 
to hold, only with a function h which depends on 5. Naturally, h{po,p 2 ] 5 ) is then given by 

Since to our knowledge the explicit form of the solution function h is not available in the 
literature, let us present some details for its derivation. For the moment, let us return to 
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standard arithmetic for and consider pj G M, z = 0,1, 2, and take Ps = Po + Pi ~ P 2 - Next 
change variables from pi to 


u = 


P2 -Po Pi- P 3 


V = 


P2 +P0 


W = 


Pi +P 3 


2 2 2 

Then 

Pi = W + U , P3 = W — U , P 2 = V + U , Po = V — U . 

Therefore, the energy constraint ^4 = 0 is equivalent to 

g{w,u) = g{v,u ), 

where, for simplicity, we set a; = 1 and then dehne 

g{v, u) = ijl — 25 cos('(; + u) — ^/l — 25 cos(u — u) = uj 2 — ojq ■ 
and thus g{w, u) = ui — UJ 3 . 


(119) 

( 120 ) 

( 121 ) 

( 122 ) 


Since g{v,TTn) = 0 for all v and n G Z, (121) is solved by any i;,u) G M if u G vrZ. This 
corresponds to the trivial solution p 2 = po mod 27r. 

Assume thus u 0 vrZ. To fix a sign convention, let us next suppose that w,v € J = (—vr, vr] 
which can always be achieved without changing u by shifting p£ by a suitably chosen integer 
multiple of 27r. Since 


, , ^cos(v — u) — cos(v + u) .^siniisinu 

g{v,u) = 25 -^- - = 45- 


W2 + UJQ 


UJ2 + Wo 


(123) 


where sin it 7 ^ 0 and oj£ > 0, (121) can only hold if sinu and sin it; have the same sign. For 
w,v £ J this is equivalent to requiring sign(ic) = sign(i;). 


To solve (121), we first note that 

g{v, it)^ = 2(1 — 2(5 cos it cos v — u}oU} 2 ), 
g{w, li)^ = 2(1 — 2(5 cos it cos w — wiwa) . 


Hence, if (121) is true, we need to have 

UJ1UJ3 = —2(5 cos it(cos w — cos v) + a;oa;2 . 

We square both sides one more time, and use the identities 

cjf u;| = 1 — 4(5 cos it cos w + 45^ cos^ w — 45^ sin^ it, 
^ 0^2 = 1 — 4(5 cos It cos V + 45^ cos^ v — 45^ sin^ u , 


(124) 

(125) 

(126) 


(127) 

(128) 


which follow from the trigonometric relation cos(a + b) cos(a — b) = cos^ a — sin^ b. 

The result can be written in terms of the variable y = cosic — cosu, and we find that 
( 121 ) implies the equation 

y‘^45‘^ sin^ it — y45{cos ii(l — W 0 W 2 ) — 25 cos 1 ;) = 0 . 


(129) 


Since sin it 7 ^ 0, this equation has exactly two solutions for y. The solution y = 0 has 
cos 1 C = cosc which yields only the solution w = v for (121) with w,v £ J since then w and 
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V must have the same sign. But then p 2 = pi, so y = 0 corresponds to the second trivial 
solution. 

This leaves only the following candidate for a non-perturbative solution: 


y = 


cos tt(l — UJQ(jJ2) — 2(5 cos V 
(5 sin^ u 


(130) 


The denominator appears to lead to singularities at = 0 and at (5 = 0. However, both 

1 — 2(5 cos u cos V — ijJoIjJ2 4(5 cos u sin^ v 


singularities are removable. Namely, (130) implies that 


y + 2cos7; = cosu- 


(5 sin^ u 


1 — 2(5 cos U cos V + (jjQljJ2 


(131) 


Since the left hand side is equal to cos re + cosf and w^v have the same sign, we again get 
only one solution which expresses re as a function of u, v. 

In summary, the above computation yields the following expression for the non-perturba¬ 
tive solution. Choose P 0 )P 2 & J = (— tt, vr]. Then also v G J and we have pi = h{po,p 2 ] <5) mod 
27r, where the non-perturbative solution is given by 


x\ P2- Po , . 

Hpo,P 2 -,o) = — - -hsign 

+ 2(5 


P2+P0 


arccos 


— cos 


P2 +Po 


sinpo + sinp 2 


1 — (5(cospo + C 0 SP 2 ) + -\/(l — 2(5 cospo)(l — 2(5 cosp 2 ) 


sm ■ 


P 2 +P 0 


(132) 


with arccos G [0, vr] denoting the principal branch. In addition, to have the correct value of 
h at the apparent discontinuity p 2 = —po, we also define sign(O) = 1: this yields h{po,p 2 ) —>■ 
TT - Po = h{po,-po), as p 2 \ -po, and h{po,P 2 ) h{Po,-Po) - 27r, as p 2 Z' -po- We 
skip the rest of the details, namely proving that the above arccos is always well-defined (i.e., 
its argument lies in [—1,1]) and that the result always provides a solution to the original 
constraint problem. 

In general, for a fixed p 2 the function h is continuous, one-to-one, and satisfies 


Lo{po) +u{h{po,P 2 )) - w(p 2 ) - w(po + Hpo,P 2 ) - P 2 ) = 0 ■ (133) 

In Fig. we display a few non-perturbative solutions, at (5 = 0.4 and at (5 = 0.5, for three 
different values of p 2 - For small 6 one finds 

Hpo,P 2 'Z) = TT-po - (5(sinpo + sinp 2 ) + 0(5 "^), (134) 

which reasonably well approximates the left hand side of Fig. The existence of the above 
limiting solution was first noted in [18] and it served as an important motivation for looking 
closer at the predictions from one-dimensional kinetic theory of phonons. 

To resolve the collisional constraints, we integrate over pi, as in the FPU case. Since 
dpZ{Po^PijP 2 ) = Z{pi) —Z{po + Pi — P 2 ) and uj'{p) = 6sinp/uj{p), this yields the collision 
operator 

C[W]{po) 

X (IU 1 IU 2 IU 3 + W0W2W3 - W0W1W3 - IU 0 IU 1 IU 2 ) , (135) 

where pi = h{po,P 2 ', S) and P 3 = Pi + Po - P 2 - 


167r(5 




dp2 


1 


U}ou}2\oo3 sin Pi — oji sinpsj 
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Figure 1: The non-perturbative solution ki = i^h{2'iTkQ,2'Kk2] S) as a function of k^ for fixed 
^2 and 8 . On the left, 8 = 0.4, and on the right, 8 = 0.5. For each <5, three solutions are 
plotted, corresponding to k 2 = 0.02 (black), 0.25 (dark grey), and 0.5 (light grey). We plot 
T X T in the extended zone scheme, the dashed lines are the boundaries of a unit cell. For 
8 < 0.5 the non-perturbative solution is smooth, while for 8 = 0.5 there are cusp singularities. 
(Reprinted from |16j l 


3 Energy transport in the kinetic theory of phonons 

3.1 Entropy and H-theorem of the phonon Boltzmann equations 

For ergodic systems, time averages will converge to ensemble averages. We have been working 
under the assumption that energy is the only ergodic variable for our particle chains. In 
addition, for infinite volume systems the canonical and microcanonical ensembles should agree, 
so to study the convergence towards equilibrium, one could inspect the evolution of the 
standard entropy density. Since it is maximized for the given energy exactly by the canonical 
ensemble, this should yield an increasing function, at least asymptotically for large times. 

Consider thus a probability measure p{q,p)ii^q(i^p for the state of the finite periodic 
system with N particles. The corresponding entropy is defined by 

'S'W = - J d^qd^ppiq,p)lnp{q,p). (136) 

For small couplings, the canonical Gibbs measure is approximately Gaussian. In general, for 
measures which are nearly Gaussian in the random vector X, and assuming that X has a mean 
p. and a covariance matrix C, we have In p{q,p) ~ ^ (ln[det(27rC')] -|- (X — p)'^C~^{X — ^)). 

This implies that the entropy of such nearly Gaussian measures satisfies 

5[p] PS ^ ln[det(27rG)] + ^ Tr 1. (137) 

Thus the entropy per particle, is then approximately equal to a constant plus 

2 ^1n(detC'). For a spatially homogeneous particle chain, one has C = where X 
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denotes the Fourier-transform. Therefore, for spatially homogeneous nearly Gaussian states 
and large enough N 


N « constant + - dklndet F{k ), 


(138) 


where T(A:) is a 2 x 2 matrix defined by 


r(l,\ _ f aU [4Qik')*,Q{k)] K[Q{k')*,p{k)]\ 


(139) 


This can be expressed in terms of the phonon fields, and hence their Wigner function VF(fc), 
using the relations given in ([^: 


= 1 (^{k)-\W{k) + W{-k)) -\{W{k) - W{-k)) 

2 V i{W{k)-W{-k)) uj{k){W{k) + W{-k)) 


(140) 


Therefore, detT(fc) = W{k)W{—k), and thus we obtain the following definition of entropy 
density for phonon systems in a kinetic scaling limit 


5 [tT] = dk\nW{k) 


(141) 


After the above preliminaries, it is satisfying to find out that the above reasoning leads 
to an entropy functional which satisfies an H-theorem for phonon Boltzmann equations. This 
holds quite generally but let us check it explicitly only for the homogeneous Boltzmann 


equations dtW{k) = C\W{k)\ where the collision operator C is given either by (107) or by 


(108). In both cases, there is a positive function G such that 


C[W]{ko)= / dk6{ko + ki-k2-k‘i) 5 {iOo + uji-(^2 -iaj^)W[G{(^()Wi] 

to 

X (iTg-^ + ITf ^ - Wg-^) . 


(142) 


Since the integrand, including both constraints, is symmetric under 0 -f-)- 1 and 2 -H- 3, and 
antisymmetric under (0,1) -H- (2, 3), we find that 


dtS[Wt] = j dkoWo^dtWtiko) 

= - / 

4 7t{0,1,2,3} 

3 


dkS{ko + ki- k2- ks)6{iJo + uji - L02 - 0J3) 


X Yl [G(Loe)Wi] + kFf ^ - W2~^ - > 0 


£=0 


(143) 


Hence, the entropy ^[Wt] is increasing along solutions of the Boltzmann equation, i.e., it 
satisfies an H-theorem. The entropy production functional satisfying dtS\Wt] = Tl[IFt], 


also has an explicit form which can be read off from the right hand side of (143). 
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3.2 Steady states 


The H-theorem (143) allows to classify all steady states of the phonon Boltzmann equation. 


Namely, suppose that W is a steady state, i.e., Wt = W is a solution to the phonon Boltzmann 
equation. Since then the left hand side of (143) is zero, the integrand defining the entropy 


production has to vanish almost everywhere on the solution manifold. 

We are only interested in nondegenerate steady states which correspond to functions W for 
which W > 0 almost everywhere. Then it is possible to apply the previous explicit solutions 
Pi = h of the collisional constraints to simplify the above problem. Namely, for any steady 
state W its inverse f{p) = W{p)~^ then necessarily satisfies 


/(Po) + f{h{po,P 2 -,S)) - f{p 2 ) - f{po-P 2 + h{po,p 2 -, 6 )) = 0, 


(144) 


for almost every po and p 2 - 

Due to the collisional constraints, any linear combination f{k) = f3u!{k) + a is obviously 


a solution to the functional equation (144). As shown in mm, these are quite generally the 
only solutions for two- and higher-dimensional crystals. 

The one-dimensional case is more intricate. For instance, h{po,p 2 ', <5) —)• n—po when <5 —)• 0. 
Hence, in the limit, any function / satisfying the symmetry requirement f{p) = —/(vr — p) is 


a solution to (144). However, we do not expect that any new solutions would appear for the 


nearest neighbour dispersion relations, i.e., if <5 > 0. This has even been rigorously proven to 
be the case for the unpinned dispersion relation with 5 = see Section 5 in HZ! for details. 


In summary, we expect that the only steady states for which W is integrable, are given 


by 




Piui{k) - p) ’ 


(145) 


where /3 > 0 and p < mincj. The last condition is necessary, since if there were a point k^ for 
which ijj{ko) = p, this would lead to a nonintegrable singularity at ko for W (the singularity 
is not integrable around x = 0 in one dimension). 


Using (142) it is straightforward to check that indeed = 0 for the functions defined 

in (145). Hence, all of them are true steady states of the phonon Boltzmann equation. The 


conservation of phonon number is reflected in the appearance of the second parameter p for 
the steady states. It is expected to be a spurious chemical potential, and one expects that 
/r —)■ 0 eventually as t —)■ oo. However, this process is not captured by the standard kinetic 
approximation which we have applied here, and resolving the issue would require a different 
approach. 

Let us also mention in passing that there are additional solutions to C[W] = 0 if one 
allows VF to be a distribution. For instance, W{k) = 5{k) formally provides a solution to 
the FPU collision operator. This is not surprising since the additional constraints given by 
such distributions can restrict the solution manifold further. In general, one should treat 
such distributional solutions with great care since it is far from obvious that the oscillatory 
terms, which are neglected in the derivation of the phonon Boltzmann equation, remain lower 
order contributions for such non-chaotic initial data. However, there are known examples 
where distributional solutions seem to play a physical role, such as in the kinetic theory of 
Bose-Einstein condensation: see m and the references therein for more details. 
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3.3 Green—Kubo formula and the linearized Boltzmann equation 


Let us now return to the Green-Kubo formula for thermal conductivity. We consider here the 
kinetic prediction for the leading harmonic part of the current-current correlation function, 
as derived in Section |1.5[ We choose the basic one-parameter canonical Gibbs measure at 
temperature T > 0 as the initial st eady state. Its limiting covariance as A —?■ 0 then has a 
Wigner function as given in (145) with /? = 1/T and /r = 0. (As discussed earlier, 

we expect the true steady states to have /r = 0 even though the steady states of the kinetic 
model can have n ^ 0.) 

In fact, the results in Sec. 1.5 imply that the evolution of the harmonic part of the 
Green-Kubo correlator is determined by the phonon Boltzmann equation linearized around 
its stationary solution /?) = l/{/3uj{k)). We begin with the following approximation 

for the Green-Kubo correlator obtained from ( |69[ ) and 

C(t;/3) « f dkMk)deWi^\k)\e=o, (146) 

Jt 

where we have written the covariance in terms of the Wigner function (k) and defined 

(po{k) = sin(27rfc). (147) 

The Wigner function is computed using expectation over the stochastic process which 
starts from a perturbation of the stationary state whose initial Wigner function converges to 
WQ^\k) = l/{/3u]{k) — e(j)o{k)) as A —>■ 0. 

Since e\(j)o{k)\ < 2TTe\k\, the Wigner function is positive for sufficiently small e > 0, 
even in the FPU-models for which w(0) = 0: this follows from the observation that it is 
always possible to find some cq > 0 such that uj{k) > co\k\ in a neighbourhood of 0. Hence, 
ht{k) = dsWl;^\k)\s=o is well defined and its initial value is given by /iq = Its 

time-evolution can be determined using the equality dtht{k) = dedtWl:‘^\k)\^=Q and the kinetic 
theory approximation ~ C[W^^^] where C denotes the appropriate phonon Boltzmann 


collision operator. The explicit form is straightforward to compute using (142) and the fact 
that 5 {ijJq+u}i—oj 2 —oj^) = 0, for any of the functions W = 

in (145). We obtain the linearized Boltzmann equation dtht = —Cht where C is the operator 
for which Ch = LW~^h with W~^ denoting multiplication by the function f3u:{k) and 


{Lf){ko) = / dk 6 {ko + ki - k 2 - k3)5{oJo + uJi - UJ 2 - CJ 3 ) 

Jjs 

3 


xYl[G{ue)We] (/o +/i -/2 -/s) , 


(148) 


£=0 


with Wi = l/{l3uj{ki)) and fe = /{ki). 

The operator L is self-adjoint, even positive, on the Hilbert space L^(T): if / G L^(T) is 
any sufficiently regular function, then by the symmetry properties of the integrand, the scalar 
product between / and Lf satisfies 


(/, L/) = ^ / dk 6 {ko + ki - k 2 - k 3 ) 6 {uJo + uii - 002 - u} 3 ) 

4 Jj{0,l,2,3} 


X II [G{ue)We] x |/o + /i - /s - /aP > 0. 


(149) 


e=o 
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In fact, the positivity of the operator L is closely related to the maximization of entropy at 
the steady state. Namely, by differentiating the entropy production (143) twice, we find 


dtdlS[wl :%=0 = 2(/, Lf ), with / = W-^ht , (150) 


where we have used the above mentioned property that the integrand vanishes whenever 
either of the two factors involving differences of W~^ is left undifferentiated. 

Since L is positive operator, any operator of the form A^LA, where A is a bounded 
operator and A^ denotes its adjoint, is also positive. Therefore, the operator L = W~^LW~^ 
is positive on L^(T), and thus the solution of the evolution equation dtht = —Lht can be 
written as ht = e~^^ho where each e“*^ is a contraction operator on L^(T). 

These properties indicate that it is more natural to study the perturbations in terms of 
ht instead of ht: any solution to dtht = —Cht provides a solution to dtht = —Lht by setting 
ht = W~^ht and vice versa. Therefore, the Boltzmann equation linearized around a steady 
state W has a solution 


ht = We-^^W-^ho, 


(151) 


for every initial perturbation ho for which W~^ho G Lp'iT). 

As mentioned above, the Green-Kubo correlation function concerns the case with /lo(^) = 
lT(/c)^ sin(27rA:) and W{k) = Thus W~^ho is equal to 


hoik)=/3 


_i sin(27rA:) 
oj{k) 


(152) 


which is a bounded function on the torus; for FPU-like models with uj{k) = 0{\k\), the 
function ho is not continuous at A: = 0 but its left and right limits exist and are finite. Hence, 
we can use the solution in (151) and conclude that the Green-Kubo correlation function at 
the steady state with temperature /3~^ has a kinetic theory approximation 


Cit- /3) « ht) = jdk Mk)W{k) ((T^'^ho^ (k) = uJ^5^{ho, e-^^ho ). (153) 

The right hand side is equal to the L^-norm \\uj'^5e~2^^ho\\'^, in particular, it is always 
positive. The question about the prediction of kinetic theory for the thermal conductivity 
at a certain equilibrium state hence boils down to the decay of the above norm under the 
semigroup e“*^. As we will show next, the two Boltzmann equations discussed above prove 
that both integrable and non-integrable decay can occur in the kinetic theory of phonons. 


3.4 Kinetic theory prediction for thermal conductivity in chains with an- 
harmonic pinning 


Suppose the positive operator L has a spectral gap of size (io > 0 above its zero eigenvalue and 
ho is orthogonal to its eigenspace of zero. Then spectral theory implies a bound \\e~^^^hoW^ < 
which is integrable over t. Hence, in this case the kinetic theory prediction is 
always a finite conductivity and its leading behaviour in A 4 can be computed using the kinetic 


approximation in (59). This yields 

poo 

^ / dr(/io,e 


—rL 


ho) = lim 

£^ 0 + , 


dr(^o,e 


= lim /3'^uj‘^6‘^{ho,{e + L) 

£-S>0 + 


-1 


ho) 


(154) 
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where the middle equality is a consequence of dominated convergence theorem. 

We have added the regulator e > 0 here to make the operator L invertible on the whole 
Hilbert space. Then the scalar product (ho, (e + L)~^ho) is finite even if ho is not orthogonal 
to the zero subspace of L. This has the benefit that one can use the full operator L, instead 
of its restriction to the orthocomplement of the zero eigenspace, and study the expectation 
of the resolvent, 


R{e) = (ho,ie + L) ^ho), e>0, (155) 

instead of the asymptotic behaviour of the semigroup. 

The above result implies that, should lim£_j.o+ h?(e) = oo, then the kinetic prediction for 
the thermal conductivity at the corresponding steady state is also inhnite (note that both 
L and ho depend on the choice of the steady state). In particular, this happens if ho is not 
orthogonal to the zero eigenspace of L: otherwise, we have R{s) > ||Po^o|P/^ oo, as e —>■ 0, 
where Po denotes the orthogonal projection to the zero eigenspace. Hence, the orthogonality 
condition mentioned in the beginning of this subsection is necessary for a finite prediction 
from kinetic theory. 

Let il{da) denote the spectral measure obtained from the spectral decomposition of L 
with respect to the vector ho- Then /2 is a positive Borel measure on the real axis whose 
support lies on the spectrum of L, on (t{L) C M+, and it is characterized by the condition 

{ho, f{L)ho) = j Jl{da)f{a) , (156) 


which holds, in particnlar, for every continuous real function /. Thus ||ho|P = f M(da), and 

1 


R(e) = 


/2(da) 
t(L) £ + a 


(157) 


If ho = 0, we have R(£) = 0 for all e, and thus in this case the kinetic theory prediction 
would be zero conductivity. If ho 7 ^ 0, such as for the above discussed particle chains, we can 
normalize the measure /I into a probability measure by dividing it by ||ho|P- In addition, the 
second derivative of the function a 1 —)• (e+a)“^ is positive, and thus it is convex. Therefore, 
we can apply Jensen’s inequality in ( 157|) to conclude that 


R(£) = 


/2(da) 1 


> 


{ho, ho)" 


(158) 


The right hand side converges to {ho,ho)‘^/{ho,Lho) as e 
ho = W(j)o, we find that 

„2-A,2{W(Po,Wct>o)^ 


(Z)||hoPe + a e{ho,ho) + {ho, Lho) 

0. Since L = W-^LW-^ and 


(159) 


forms a lower bound for the kinetic theory prediction of thermal conductivity in the particle 
chains. _ 

By (1491, V’ can be a zero eigenvector of L if and only if it belongs to L^(T) and 


f = W satisfies (144), i.e., it is a collisional invariant. We have argued (and even 
proven for 5 = 1/2) that all collisional invariants are linear combinations of cn and 1. This 
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implies that the zero subspace of L is spanned by 1 and W. For the particle chains, 
ho{k) = W{k)4)Q{k) = u'{k)/{2 t: 5up'P) where uj' denotes the derivative of uj. Therefore, 
by periodicity, jjdikho{k){aW{k) + 6) = 0 for all o, 6 G C, and thus ho is indeed orthogonal 
to the proposed zero subspace of L. 

Let us also sketch a possible proof for the spectral gap of L. Cons ider an arbitrary p and 
denote / = W~^p. We begin from the integral representation (149) for {f,Lf) = {ip, Lip). 
Let us mollify the singularities in the integrand by choosing a suitable regularizing function 
—>■ [0,1] and then using the bound 1 > ^{k) inside the integrand. The purpose of 
4> is to regularize the collision cross section of the model, so that in the remaining integral 
one can use the explicit solution and define an operator L' such that (/, Lf) > (/, L'f) and 
L' = V — K where IL is a positive multiplication operator and iL is a self-adjoint integral 
operator. If we can then tune so that V and XjY are both bounded and B = V~^PKV~^P 
is a compact operator, it follows that 1 — B has a spectral gap and its zero subspace consists 
of those functions (p for which V~^P(p is a collisional invariant. Then for every / which is 
orthogonal to the collisional invariants, we have (/, L/) > {f,L'f) > (j^||/|P for some 5' > 0 
independent of /. This estimate would then prove a gap also for L in the pinned case: then 
W and \/W are both bounded functions and there is (5o > 0 such that {ip, Lip) > (joUV'lP for 
every ip for which W~^ip is orthogonal to the collisional invariants. 

The above procedure is a variant of the standard argument used to prove a gap in kinetic 
theory when the “relaxation time” determined by the multiplication operator V is bounded 
from above and below (see the next section for discussion about the relaxation time ap¬ 
proximation). The problem about using the standard argument directly for the linearized 
operator of the pinned chains lies in the non-integrable singularity related to the Jacobian of 
the change of variables which resolves the energy constraint. This singularity is the same one 
which makes the total collision cross section—obtained by neglecting all terms which contain 


Wo in (135)—infinite: since in the integrand ps = pi — y when p 2 = po + y, the integrand 
has a nonintegrable singularity of strength at least at p 2 = Po- Thus both the total 

collision cross section and the relaxation time function V are formally infinite for all ko- The 
sign however is positive, and thus by reducing the strength of the collisions one can aim at 
approximating the linearized operator from below by an operator of the “standard form” with 
a finite relaxation time function. However, let us not try to complete the argument here but 
instead focus on its implications on the thermal conductivity. 

To summarize, the above argument strongly indicates that the kinetic theory prediction 
for the thermal conductivity at temperature for chains with anharmonic pinning is a 
finite non-zero value, and for small A 4 we should have 


Kirp- W^^S^{ho,L-^ho), 


(160) 


where “L~^ho” denotes the unique (p G L^(T'^) for which L(p = ho (such (p can be found if L 
has a spectral gap and ho is orthogonal to the zero eigenspace of L, as we have argued above). 
Inserting W{k) = {j5uj{k))~^ and using the definitions of L and ho this can be simplified to 

k(/3“^) p~‘^uj'^6^{uj~^(po, L~^{uj~‘^(po)) , (161) 
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where the operator L is explicitly given by 


{L'if;){ko) = I dk 6 {ko + ki - k 2 - k 3 ) 6 {u}o + oji - lo 2 - 003 ) 

4 Jt3 


X + i’l - 4^2 - 'ips) ■ 


e=o 


Thus the dependence on the temperature and coupling factorizes: 

k (/ 3 -^) l3^Xfu^C{6) , 

where the constant C{5) is a function of the harmonic pinning parameter 5 only, 


,8 


C{S) = 6 '^-{i' (Po,Lq {u (po)), with 


{Loi>)iko) = 


dp 2 

. 2-71 


sinpi smp 3 


121 


1^3 


-1 3 


n ^ (V'O + V’l - ^^2 - -03) 


^0 
i=0 ^ 


( 162 ) 

(163) 

(164) 

(165) 


wher e -ipi = • 0 (P£/( 2 vr)), po = 2 -Kko, ps = pi + pQ - p 2 , and pi = h{po,P 2 ]S), as defined in 
(132). In addition, we have used here 12 to denote the normalized dispersion relation with 


w = 1, that is, in the above formulae z^(/c) = -^1 — 25cos(27r/c) and = u{pi/{27t)). The form 
is amenable for numerical inversion of the operator Lq, by choosing a suitable orthonormal 
basis for the subspace of L^(T) which consists of vectors orthogonal to 1 and ui. 


The Jensen inequality lower bound given in (159) implies that 




(166) 


9 {4)q,Lo(Pq) 

Using the symmetrized form in (149) for {(j)Q, LQ(j)o) and changing the integration variable 
from ko to po = 27rA;o then yields 
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S-^Ci6) > - 


dpo 


■ 2 \ 2 

sm Po 


dpo / dp 2 


sinpi smp 3 


I 2 l 


V 3 


-1 3 


n ^ IV’O + V’l - "02 - '031' 




-1 


(167) 


where 0^ = sinp^. When J —>■ 0, we have —>■ 1, pi —)• tt — po, and pi —)• vr — p 2 in the 
above. The remaining integrals can be computed explicitly, and the limit of the right hand 
side is found to be equal to 7r^/36 ~ 0.274. The numerical inversion of the full operator in 
the J —>• 0 limit in |16j resulted in the value 0.2756 which is very close to the above Jensen 


bound. In addition, evaluation of the right hand side of (167) by numerical integration shows 


that it depends only weakly on J, decreasing to 0.2 at J = 0.3. However, the bound becomes 
ineffective for larger values of J, going to zero as <5 — 

Therefore, the formula 0.274 provides a fairly good approximation for the lower 


bound in (159) for the kinetic theory prediction of the thermal conductivity, at least for small 


enough 6 . This approximation was compared in [T^ to the thermal conductivity measured in 
numerical simulations of large finite chains with boundary thermostats. The numerical simu¬ 
lations where performed with tJ = and A 3 = 0, and instead of A 4 —?■ 0 with fixed f5 one 
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Figure 2: Measured T‘^k(T) compared against the kinetic theory computation for small 5 and 
T (straight line). The data points are for T = 0.1 (□), 0.4 (o), and 4(A). (Reprinted from 

m) 


considers /3 —>■ oo with fixed A 4 . The two limits can be connected by a straightforward scaling 
argument which allows to compare the kinetic prediction with the conductivity observed in 
the simulations (see m Section 2] for details). As mentioned above, the Jensen inequality 


lower bound given in (159) is a good approximation of the numerically inverted value for small 
5. In Fig. we have given the comparison between the values from the numerical particle 
simulations and the Jensen bound. The agreement is surprisingly good and seems to indicate 
that for this model the kinetic theory gives good description of the dominant effects affecting 
thermal conduction in the pinned anharmonic chains. 


3.5 Anomalous energy conduction in the kinetic theory of FPU chains 

The kinetic approximation to the Green-Kubo correlation function in the FPU chains has 
been studied rigorously in m- As mentioned above, the FPU chains have 5 = and in 
m the frequency normalization was chosen as 11^/2 since then oj{k) = |sin(7r/c)| has no 
prefactors. For easy comparison, let us make this choice also in this subsection for the FPU 
chains: if needed, the scale u can be reintroduced in the results just as was done in the 
previous subsection. 


With this choice and using (148), the Boltzmann collision operator linearized around a 
steady state W{k) = /3“^a;(A:)“^ becomes 


{Lf){ko) = 97rXl{2/j3r / dk S{ko + ki - k 2 - k3)d{uJo + ui - UJ 2 - UJ 3 ) 

X (/o + /l - /2 - /s) , 


(168) 
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and it is related to the kinetic theory prediction of the Green-Kubo correlation function by 

Cit; P) ^ ^{ho, e~^^ho) , (169) 

where L = /3'^ujLuj and ho{k) = 2 cos(7rA:)//3 for k G [0,1). Choosing a normalization somewhat 
different from the previous subsection, we may rewrite this in a dimensionless form by defining 
^^{k) = i cos(7r/c), 0 < A: < 1, and Lq = ^(/3/2)^t<;La;. This yields 

C{t; 13) « , (170) 

where the constant c = 12^A|/3“^7r“^ and Lq = loLqlo. In m, this was called the kinetic 
conjecture and the above form coincides with the one given in m Eq. (1.18)] after a change 
of variables from k G [0,1) to p G [0, 27r). 

Explicitly, the operator Lq can be written for p G I = [0, 27r) as 

(Lo/)(po/(27r)) = / dp5(po+pi -/C 2 - A:3)5(a;o+wi-a ;2 - cus) 

X (/o + /i - /2 - h) . (171) 

Now the constraints can be explicitly integrated using the results mentioned in Sec. 
yielding 


where 


Lq = V — A , with 

{V'tl>){k) = V{27rk)‘ip{k) and {A'i(^){k) = f dk'27rK{2TTk,27rk')'tp{k'), 

Jo 


V{p) = / dp'K 2 {p,p '), K{p,p') = 2K2 {p,p') - Ki{p,p') , with 

= and K,(p.p-)= " 


V^F_(p,p') 


\/F+(p,p') 


(172) 

(173) 

(174) 

(175) 


In this formula, the multiplication operator arises from the ^^fQ^’-term and the Ki term in 
the integral kernel K from the “/i”-term. By symmetry, the contributions from the “/2”- 
and “/a” -terms are equal, each contributing one Er2-term to the integral kernel. It probably 
comes as no surprise that the precise analysis of the semigroup generated by Lq gets rather 
technical. However, this has been done in m, and let us only repeat the main conclusions 
from the analysis here. 

Unlike for the onsite anharmonic perturbation, the relaxation time for the EPU models is 
finite, and as the above formula shows, the operator Lq can be written in the standard form 
V — A where V = oj'^V is a multiplication operator and A = ujAoj is an integral operator. It 
is proven in m Lemma 4.1] that the function V{k) is continuous and can be bounded from 
above and below by j sin(7rfc)j^/^. In particular U(0) = 0 and, consequently, the operator Lq 
has no spectral gap. 

In kinetic theory, it is a common practice to use the relaxation time approximation to 
approximate the linearized Boltzmann evolution. In our case, this amounts to dropping the 
operator A, i.e., approximating 

{tp, ijj) « • (176) 
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Since V'(O) = 1) the decay of the relaxation time approximation is now entirely determined 
by the values of the potential near zero. The above bounds imply that V{k) = 0{\k\^/‘^), 
and thus the relaxation time approximation predicts C{t; /3) = for large t. (To our 

knowledge, this decay of the relaxation time approximation was first as derived in m-) 

However, the contribution arising from adding the integral operator A is also singular, and 
more careful study is required to conclude that the relaxation time prediction continues to 
hold for the full semigroup. Fortunately, the above straightforward estimate gives the correct 
decay speed: it was shown in m that the resolvent of the semigroup satisfies 
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(177) 


where the first term is identical to the relaxation time approximation, and behaves as 
for small e > 0. The second and third term are for any e > 0. Although also 

this second contribution is divergent, the first term is dominant, and thus we confirm the 
prediction of the relaxation time approximation in this particular case. 

Even the exact asymptotics can be found for this particular choice. Applying m Corollary 
2 . 6 ] we find that the kinetic theory predicts for times with t(A 4 // 3 )^ S> 1 , and for sufficiently 
small couplings As and A 4 , 

C{t;P)^CoP-Hxlt)-l . (178) 


Here Cq = co('7rl2“^)^/®/(27rr(2/5)) is an explicit numerical constant. Evaluation of the 
Gamma-function and the integrals defining the constant “cq” in (6.14) and (4.7) of [IT] yields 
a numerical approximation Cq ~ 0.00386. 

Of course, for much longer than kinetic times, the terms neglected in the derivation 
of the Boltzmann equation might become important and alter the asymptotic decay. The 
above results also imply that, on the kinetic time scale, the energy spread is superdiffusive: 
the quadratic energy spread observable discussed in Sec. 

S{t) = 

The energy spread has been analysed in more detail in |22j . The authors study the time 
evolution of inhomogeneous perturbations around a given thermal equilibrium state. This 
results in an evolution equation which corresponds to the phonon Boltzmann evolution where 
the nonlinear collision operator has been replaced by the above linearized operator. Explicitly, 
the perturbation f{t,x,k), defined via W = {1 + f)/{f3uj), evolves then by 


1.5 should then be increasing as 


dtf{t, X, k) + ^J{k)dxf{t, X, k) = {Lf{t, x, ■)){k). 


(179) 


The result concerns L^-integrable initial data varying at a scale with e small. It is shown 
that for sufficiently long times, the solution then first thermalizes in the /c-variable, which 
by the definition of / implies that it becomes independent of k. Diffusive relaxation in the 
spatial variable would then mean that at a time scale the perturbation follows the heat 
equation dtf + k(—A)/ = 0, with k > 0. However, this does not occur here: instead, it is 
shown that at the time scale the perturbation satisfies a fractional diffusion equation 

dtf + «:(—A)^/^/ = 0 , with k > 0 . 

This corresponds to a superdiffusive relaxation of the initial perturbation. Moreover, the 
fractional diffusion spreads local perturbations at time t only up to distances This is 
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in apparent contradiction with the earlier claim that S{t) = which would indicate 

that the spatial spread occurs at a speed i.e., faster than predicted by the fractional 

diffusion equation. 

The resolution lies in the tail behaviour of solutions to the fractional diffusion equation. 
Let us conclude the section with a somewhat heuristic argument which would explain the 
above results. By using Fourier-transform and a simple scaling argument, one finds that the 
solution to dtf+K{—A)^/^f = 0, with an initial data /o, is given by a convolution of /o with an 
integral kernel Kt which satisfies a scaling relation Kt{x) = t~^F{x/t'^) where p = 5/8. Unlike 
the Gaussian heat kernel in ( |46[ ), the function F is merely polynomially decreasing, with F{y) 
decaying as for large \y\ (the power 13/5 is obtained by evaluating 1/p + 1). Thus 

Jdyy'^F{y) = oo, and the quadratic spread from fractional diffusion becomes immediately 
infinite, even if it is initially finite. 

However, on the microscopic scale, the velocities of the ballistic harmonic evolution are 
bounded from above, and this will eventually cut off the above decay, and change it from 
the above powerlaw decay to an exponentially fast decay at spatial distances 0 {t) from the 
source. Therefore, we would expect that the true microscopic distribution of an initially local 
perturbation at a large time t is 0 {t~P) for distances |x| = 0 {tP), it is for |x| 

between O^F) and 0(t), and it becomes exponentially decreasing for distances larger than 
0{t). For such functions, the value of S{t) is entirely dominated by the midscale powerlaw 
tail which yields a term just as we obtained from the kinetic prediction for the 

Green-Kubo correlation function using the linearized Boltzmann equation. 

4 Concluding remarks 

As the last two explicit examples show, kinetic theory of phonons is capable of uncovering 
detailed information about the decay of time-correlations and, via the Green-Kubo formula, 
about the thermal conductivities of classical particle chains. This is somewhat surprising, 
considering the various mathematical problems and uncertainties discovered along the way to 
the phonon Boltzmann equation and even in its analysis. However, the agreement between 
the kinetic prediction and numerical simulations for thermal conductivities in chains with 
anharmonic pinning, and the discovery of anomalous energy transport by fractional Brownian 
motion from the linearized Boltzmann equation in the FPU-/3 chains, present a strong case 
in favour of looking for further applications of the phonon Boltzmann equations, even in the 
somewhat degenerate one-dimensional case. After all, it is at present one of the very few 
general tools which allow computing the dependence of the thermal conduction properties 
on the parameters of the microscopic evolution directly, without introduction of additional 
fitting parameters. 

We have also seen that some care is needed in the application of the phonon Boltzmann 
equation. Most importantly, the equation is closely tied to the scales on which the free 
streaming of phonons begins to alter its character due to the collisions. Hence, it describes 
the evolution up to the kinetic time-scale only, and it is possible that further changes are 
found at larger time-scales. Nevertheless, since the H-theorem implies that solutions to the 
kinetic equation push the system towards thermal equilibrium states, drastic changes in the 
character of the evolution should be the exception, not the rule. 

The evaluation of the decay of Green-Kubo correlation functions is one of the robust 
applications of kinetic theory, and the thermal conductivity obtained by the perturbation 
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procedure recalled here should in general yield its leading behaviour in the limit of weak 
perturbations. However, as a warning about the standard procedure, let us stress that the 
commonly used relaxation time approximation of the linearized collision operator is only an 
order-of-magnitude estimate of the real kinetic prediction which involves inverting the full lin¬ 
earized operator. In the pinned case, we found that the relaxation time approximation predicts 
zero thermal conductivity. This turned out to be misleading since already the straightforward 
Jensen bound for the full inverse proves that the kinetic prediction for the conductivity is 
non-zero. In contrast, for the anomalous conduction in the FPU-/3 chain the relaxation time 
approximation does capture the correct asymptotic decay of the kinetic prediction. 

One major open question in the kinetic theory of phonons, and in fact of nonlinearity 
perturbed systems in general, is the precise manner of handling spatially inhomogeneous 
perturbations. The standard Boltzmann transport term, appearing on the left hand side 
of (179), might require adjusting to capture all effects relevant to the transport. One such 
example often found in the literature is an addition of a Vlasov-type term. For which systems, 
under which time-scales, and for which initial data such corrections are necessary, remains 
unresolved at the moment. 

One benefit from a better understanding of the behaviour of inhomogeneous perturbations 
could be a first-principles derivation of fluctuating hydrodynamics for these systems, including 
the precise dependence of its parameters on the microscopic evolution, in the limit of weak 
couplings. The application of fluctuating hydrodynamics to the transport in one-dimensional 
particle chains has been discussed in [23] and reviewed in |23|. It appears to be the first 
model which is able to describe the anomalous transport in one-dimensional particle chains 
fully in agreement with computer simulations of the spread of localized perturbations. Con¬ 
necting it directly to the microscopic dynamics would be a breakthrough in understanding 
the microscopic origin and precise nature of transport in crystalline structures, such as the 
present particle chains. Kinetic theory, the phonon Boltzmann equation in particular, could 
well provide some of the missing steps into this direction. 


Acknowledgments 

I am most grateful to Herbert Spohn for his comments and suggestions for improvements. 
Most of the discussion here is based on his works and on our joint collaborations. The related 
research has been made possible by support from the Academy of Finland, and benefited 
from the support of the project EDNHS ANR-14-CE25-0011 of the French National Research 
Agency (ANR). I am also grateful to Matteo Marcozzi and Alessia Nota for their comments 
on the manuscript. 


References 

[1] H. Spohn, The phonon Boltzmann equation, properties and link to weakly anharmonic 
lattice dynamics, J. Stat. Phys. 124(2-4) (2006) 1041-1104. 

[2] J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids. 
Oxford University Press, London, 1967. 

[3] M. Butz, Kinetic limit for Wave Propagation in a Continuous, Weakly Random Medium: 
Self-averaging and Convergence to a Linear Boltzmann Equation. PhD thesis, Technische 
Universitat Miinchen, 2015. 


45 



[4] J. Lukkarinen and H. Spohn, Kinetic limit for wave propagation in a random medium, 
Arch. Ration. Mech. Anal. 183(1) (2007) 93-162. 

[5] L. Erdos, M. Salmhofer, and H.-T. Yau, Quantum diffusion of the random Schrodinger 
evolution in the scaling limit I. The non-recollision diagrams, Acta Math. 200(2) (2008) 
211-277. 

[6] L. Erdos and H.-T. Yau, Linear Boltzmann equation as the weak coupling limit of a 
random Schrodinger equation, Commun. Pure Appl. Math. 53(6) (2000) 667-735. 

[7] T. Chen, Localization lengths and Boltzmann limit for the Anderson model at small 
disorders in dimension 3, J. Stat. Phys. 120 (2005) 279-337. 

[8] P. Gerard, P. A. Markowich, N. J. Mauser, and E. Paupaud, Homogenization limits and 
Wigner transforms, Commun. Pure Appl. Math. 50 (1997) 323-379. 

[9] A. Mielke, Macroscopic behavior of microscopic oscillations in harmonic lattices via 
Wigner-Husimi transforms. Arch. Ration. Mech. Anal. 181 (2006) 401-448. 

[10] L. Harris, J. Lukkarinen, S. Teufel, and E. Theil, Energy transport by acoustic modes of 
harmonic lattices, SIAM J. Math. Anal. 40(4) (2008) 1392-1418. 

[11] D. Benedetto, E. Castella, R. Esposito, and M. Pulvirenti, From the N-body Schrodinger 
equation to the quantum Boltzmann equation: a term-by-term convergence result in the 
weak coupling regime, Commun. Math. Phys. 277(1) (2008) 1-44. 

[12] J. Lukkarinen and H. Spohn, Weakly nonlinear Schrodinger equation with random initial 
data. Invent. Math. 183(1) (2011) 79-188. 

[13] O. E. Lanford, J. L. Lebowitz, and E. H. Lieb, Time evolution of infinite anharmonic 
systems, J. Stat. Phys. 16(6) (1977) 453-461. 

[14] P. Butta, E. Caglioti, S. Di Ruzza, and C. Marchioro, On the propagation of a perturba¬ 
tion in an anharmonic system, J. Stat. Phys. 127(2) (2007) 313-325. 

[15] J. Lukkarinen and M. Marcozzi, Wick polynomials and time-evolution of cumulants, 
ArXiv e-print (2015). arXiv.org;1503.05851. 

[16] K. Aoki, J. Lukkarinen, and H. Spohn, Energy transport in weakly anharmonic chains, 
J. Stat. Phys. 124 (2006) 1105-1129. 

[17] J. Lukkarinen and H. Spohn, Anomalous energy transport in the FPU-fd chain, Commun. 
Pure Appl. Math. 61(12) (2008) 1753-1786. 

[18] R. Lefevere and A. Schenkel, Normal heat conductivity in a strongly pinned chain of 
anharmonic oscillators, J. Stat. Mech. 2006(02) (2006) L02001. 

[19] H. Spohn, Collisional invariants for the phonon Boltzmann equation, J. Stat. Phys. 124 
(2006) 1131-1135. 

[20] X. Lu, The Boltzmann equation for Bose-Einstein particles: Regularity and condensa¬ 
tion, J. Stat. Phys. 156(3) (2014) 493-545. 


46 



[21] A. Pereverzev, Fermi-Pasta-Ulam (3 lattice: Peierls equation and anomalous heat con¬ 
ductivity, Phys. Rev. E 68(5) (2003) 056124. 

[22] A. Mellet and S. Merino-Aceituno, Anomalous Energy Transport in FPU-/3 Chain, J. 
Stat. Phys. 160(3) (2015) 583-621. 

[23] H. Spohn, Nonlinear fluctuating hydrodynamics for anharmonic chains, J. Stat. 
Phys. 154(5) (2014) 1191-1227. 

[24] H. Spohn, Fluctuating hydrodynamics approach to equilibrium time correlations for an¬ 
harmonic chains, ArXiv e-print (2015). arXiv.org:1505.05987. 


47 



